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The  main  results  obtained  and  published  during  the  period  covered  by  this  report,  August  1988 
-  July  1989,  are  described  below  together  with  references  given  to  the  corresponding 
publication. 

1.  The  Interacting  Multiple  Model  Algorithm  for  Systems  with  Markovian  Switching  Coefficients, 
(Henk  A.  Blom  and  Yaakov  Bar-Shalom,  IEEE  Transactions  on  Automatic  Control  Vol.  33, 

No.  8,  August  1988) 

An  important  problem  in  filtering  for  linear  systems  with  Markovian  switching  coefficients 
(dynamic  multiple  model  systems)  is  the  one  of  management  of  hypotheses,  which  is  necessary 
to  limit  the  computational  requirements.  A  novel  approach  to  hypotheses  merging  has  been 
developed  for  this  problem.  The  novelty  lies  in  the  timing  of  hypotheses  merging.  When 
applied  to  the  problem  of  filtering  for  a  linear  system  with  Markovian  coefficients  this  yields  an 
elegant  way  to  derive  the  interacting  multiple  model  (IMM)  algorithm.  Evaluation  of  the  IMM 
algorithm  makes  it  clear  that  it  performs  very  well  at  a  relatively  low  computational  load.  These 
results  imply  a  significant  change  in  the  state  of  the  art  of  approximate  Bayesian  filtering  for 
systems  with  Markovian  coefficients. 

2.  Failure  Detection  Via  Recursive  Estimation  for  a  Class  of  Semi-Markov  Switching  Systems, 

(L.  Campo,  P.  Mookeijee  and  Y.  Bar  Shalom,  Proceedings  1988  IEEE  CPC.  Austin,  Texas) 

An  area  of  current  interest  is  the  estimation  of  the  state  of  discrete-time  stochastic  systems  with 

parameters  which  may  switch  among  a  finite  set  of  values.  The  parameter  switching  process  of 

interest  is  modeled  by  a  class  of  semi-Markov  chains.  This  class  of  processes  is  useful  in  that 

it  pertains  to  many  areas  of  interests  such  as  the  failure  detection  problem,  the  target  tracking 

problem,  socio-economic  problems  and  in  the  problem  of  approximating  nonlinear  systems  by 

a  set  of  linearized  models.  It  is  shown  in  this  paper  how  the  transition  probabilities,  which 

govern  the  model  switching  at  each  time  step,  can  be  inferred  via  the  evaluation  of  the 

conditional  distribution  of  the  sojourn  ,,ne  Following  this,  a  recursive  state  estimation 

algorithm  for  dynamic  systems  with  no*  oservations  and  changing  structures,  which  uses  . 

•* 

the  conditional  sojourn  time  distribution,  is  derived  and  and  applied  to  a  failure  detection 
problem. 
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3.  Distributed  Adaptive  Estimation  with  Probabilistic  Data  Association,  (K.C.  Chang  and  Y. 
Bar-Shalom,  Automatica.  Vol.  25,  No.  3,  pp.  359-369,  1989) 

The  probabilistic  data  association  filter  (PDAF)  estimates  the  state  of  a  target  in  a  cluttered 
environment.  This  subopdmal  Bayesian  approach  assumes  that  the  exact  target  and 
measurement  models  are  known.  However,  in  most  practical  applications,  there  are  difficulties 
in  obtaining  an  exact  mathematical  model  of  the  physical  process.  In  this  paper,  the  problem  of 
estimating  target  states  with  uncertain  measurement  origins  and  uncertain  system  models  in  a 
distributed  manner  is  considered.  First,  a  scheme  is  described  for  local  processing,  then  the 
fusion  algorithm  which  combines  the  local  processed  results  into  a  global  one  is  derived.  The 
algorithm  can  be  applied  for  tracking  a  maneuvering  target  in  a  cluttered  and  low  detection 
environment  with  a  distributed  sensor  network. 

4.  An  Adaptive  Dual  Controller  for  a  MIMO-ARMA  System,  (P.  Mookerjee  and  Y.  Bar-Shalom, 
IEEE  Transactions  on  Automatic  Control.  Vol.  34,  No.  7,  July  1989) 

An  explicit  adaptive  dual  controller  has  been  derived  for  a  multiinput  multioutput  ARMA 
system.  The  plant  has  constant  but  unknown  parameters.  The  cautious  controller  with  a 
one-step  horizon  and  a  new  dual  controller  with  a  two-step  horizon  are  examined.  In  many 
instances,  the  myopic  cautious  controller  is  seen  to  turn  off  and  converges  very  slowly.  The 
dual  controller  modifies  the  cautious  control  design  by  numerator  and  denominator  correction 
terms  which  depend  upon  the  sensitivity  functions  of  the  expected  future  cost  and  avoids  the 
turn-off  and  slow  convergence.  Monte-Carlo  comparisons  based  on  parametric  and 
nonparametric  statistical  analysis  indicate  the  superiority  of  the  dual  controller  over  the  cautious 
controller. 
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5.  Time- Reversion  of  a  Hybrid  State  Stochastic  Difference  System,  (Henk  A.P.  Blom  and 
Yaakov  Bar-Shalom,  Proc.  1989  rEEE  Intni.  Conf.  on  Control  &  Applications.  Jerusalem, 
Israel,  April  1989  to  appear  in  IEEE  Trans.  Info.  Theory.  1990) 

This  paper  develops  the  reversion  in  time  of  a  stochastic  difference  equation  in  a  hybrid  space, 
with  a  Markovian  solution.  The  reversion  is  obtained  by  a  martingale  approach,  which 
previously  led  to  reverse  time  forms  for  stochastic  equations  with  Gauss-Markov  or  diffusion 
solutions.  The  reverse  time  equations  follow  from  a  particular  non-canonical  martingale 
decomposition,  while  the  reverse  time  equations  for  Gauss-Markov  and  diffusion  solutions 
followed  from  the  canonical  martingale  decomposition.  The  need  for  the  non-canonical 
decomposition  stems  from  the  hybrid  state  space  situation.  The  non-Gaussian  discrete  time 
situation  leads  to  reverse  time  equations  that  incorporate  a  Bayesian  estimation  step. 

6.  A  New  Controller  for  Discrete-Time  Stochastic  Systems  with  Markovian  Jump  Parameters,  (L. 
Campo  and  Y.  Bar-Shalom,  11th  IFAC  World  Congress.  Tallinn,  USSR,  Aug.  1990 

A  realistic  stochastic  control  problem  for  hybrid  systems  with  Markovian  jump  parameters  may 
have  the  switching  parameters  in  both  the  state  and  measurement  equations.  Furthermore,  both 
the  system  state  and  the  jump  states  may  not  be  perfectly  observed.  Prior  to  this  work  the  only 
existing  implementable  controller  for  this  problem  was  based  upon  a  heuristic  multiple  model 
partitioning  (MMP)  and  hypothesis  pruning.  In  this  paper  a  stochastic  control  algorithm  for 
stochastic  systems  with  Markovian  jump  parameters  was  developed.  The  control  algorithm  is 
derived  through  the  use  of  stochastic  dynamic  progamming  and  is  designed  to  be  used  for 
realistic  stochastic  control  problems,  i.e.,  with  noisy  state  obeservations.  The  state  estimation 
and  model  identification  is  done  via  the  recently  developed  Interacting  Multiple  Model 
algorithm.  Simulation  results  show  that  a  substantial  reduction  in  cost  can  be  obtained  by  this 
new  control  algorithm  over  the  MMP  scheme. 
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7.  From  Piecewise  Deterministic  To  Piecewise  Diffusion  Markov  Processes,  (Henk  A.P.  Blom, 
Proc.  IEEE  CPC  1988) 

Piecewise  Deterministic  (PD)  Markov  processes  form  a  remarkable  class  of  hybrid  state 
processes  because,  in  contrast  to  most  other  hybrid  state  processes,  they  include  a  jump 
reflecting  boundary  and  exclude  diffusion.  As  such,  they  cover  a  wide  variety  of  impulsively 
or  singularly  controlled  non-diffusion  processes.  Because  PD  processes  are  defined  in  a 
pathwise  way,  they  provide  a  framework  to  study  the  control  of  non-diffusion  processes  along 
same  lines  as  that  of  diffusions.  An  important  generalization  is  to  include  diffusion  in  PD 
processes,  but,  as  pointed  out  by  Davis,  combining  diffusion  with  a  jump  reflecting  boundary 
seems  not  possible  within  the  present  definition  of  PD  processes.  This  paper  presents  PD 
processes  as  pathwise  unique  solutions  of  an  Ito  stochastic  differential  equation  (SDE),  driven 
by  a  Poisson  random  measure.  Since  such  an  SDE  permits  the  inclusion  of  diffusion,  this 
approach  leads  to  a  large  variety  of  piecewise  diffusion  Markov  processes,  represented  by 
pathwise  unique  SDE  solutions. 
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nrerglag.  When  applied  (o  Ibe  problem  of  tio(  for  a  linear  system  with 
Markovian  coefficient!  (his  yields  an  t  .  way  to  derive  the  Interacting 
multiple  node!  (IMM)  algorithm.  E>.  atiosi  of  Ibe  IMM  algorithm 
makes  it  dear  that  it  performs  very  weil  at  a  relativeiy  low  compotational 
load.  These  results  imply  a  significant  change  in  the  state  of  tbe  art  of 
approximate  Bayesian  filtering  for  systems  with  Markovian  coefficients. 

I.  Introduction 

In  this  contribution  we  present  a  novel  approach  to  the  problem  of 
filtering  for  a  linear  system  with  Markovian  coefficients 

*,  =  a(«,)j c,.,  +  6(«,)w,  ft) 

with  observations 
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y,  =  h(6,)xl  +  g(6,)u,  (2) 

9,  is  a  finite  state  Markov  chain  taking  values  in  { 1 ,  •  •  • ,  /V}  according  to 
a  transition  probability  matrix  H,  and  w,,  u,  arc  mutually  independent 
white  Gaussian  processes.  The  exact  filter  consists  of  a  growing  number 
of  linear  Gaussian  hypotheses,  with  the  growth  being  exponential  with  the 
time.  Obviously,  for  filtering  we  need  recursive  algorithms  whose 
complexity  does  not  grow  with  time.  With  this,  the  main  problem  is  to 
avoid  the  exponential  growth  of  the  number  of  Gaussian  hypotheses  in  an 
efficient  way. 

This  hypotheses  management  problem  is  also  known  for  several  other 
filtering  situations  (10],  (5],  [6],  (9],  and  [4],  All  these  problems  have 
stimulated  during  the  last  two  decades  the  development  of  a  large  variety 
of  approximation  methods.  For  our  problem  the  majority  of  these  art 
techniques  that  reduce  the  number  of  Gaussian  hypotheses,  by  pruning 
and/or  merging  of  hypotheses.  Well-known  examples  of  this  approach  art 
the  detection  estimation  (DE)  algorithms  and  the  generalized  pseudo 
Bayes  (GPB)  algorithms.  For  overviews  and  comparisons  see  [14],  [7], 
(12),  and  [17].  None  of  the  algorithms  discussed  appeared  to  have  good 
performance  at  modest  computational  load.  Because  of  that,  other 
approaches  have  been  also  developed,  mainly  by  way  of  approximating 
the  model  (1),  (2).  Examples  are  the  modified  multiple  model  (MM) 
algorithms  [20],  [7],  the  modified  gain  extended  Kalman  (MGEK)  filter  of 
Song  and  Speyer  [13],  [7],  and  residual  based  methods  [19],  [2],  These 
algorithms,  however,  also  lack  good  performance  at  modest  computa¬ 
tional  load  in  too  many  situations.  In  view  of  this  unsatisfactory  situation 
and  the  practical  importance  of  better  solutions,  the  filtering  problem  for 
the  class  of  systems  (1),  (2)  needed  further  study. 

One  hem  that  has  not  received  much  attention  in  the  past  is  the  timing  of 
hypotheses  reduction.  It  is  common  practice  to  reduce  the  number  of 
Gaussian  hypotheses  immediately  after  a  measurement  update.  Indeed,  on 
first  sight  there  does  not  seem  to  be  a  better  moment.  However,  in  two 
recent  publications  [3],  [1],  this  point  has  been  exploited  to  develop, 
respectively,  the  so-called  IMM  (interacting  multiple  model)  and  AFMM 
t adaptive  forgetting  through  multiple  models)  algorithms.  The  latter 
exploits  pruning  to  reduce  the  number  of  hypotheses,  while  the  IMM 
exploits  merging.  The  IMM  algorithm  was  the  reason  for  a  further 
evaluation  of  the  timing  of  hypotheses  reduction.  A  novel  approach  to 
hypotheses  merging  is  presented  for  a  dynamic  MM  situation,  which  leads 
to  an  elegant  derivation  of  the  IMM  algorithm.  Next  Monte  Carlo 
simulations  are  presented  to  judge  the  state  of  the  art  in  MM  filtering  after 
the  introduction  ol  the  IMivi  algorithm. 

n.  Timing  of  Hypotheses  Reduction 

To  show  the  possibilities  of  timing  the  hypothesis  reduction,  we  start 
with  a  filter  cycle  from  one  measurement  update  up  to  and  including  the 
next  measurement  update.  For  this,  we  take  a  cycle  of  recursions  for  the 
evolution  of  the  conditional  probability  measure  of  our  hybrid  state 
Markov  process  (x,,  $,).  This  cycle  reads  as  follows: 

Mixing 

/»{«,_,[ r,.,}  — -  P{e,\r,.,)  0) 
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if  P{8, |  Y,-i}  =0  prune  hypothesis  8,, 

Muiaf 

n.i  1  -  —  p[jr, n  ,l  (4) 

Evotwooe 

pu,-,i«,.  r,-, i  - — —  />[*,{«„  y,., ]  (5) 

B*ye* 

P{8,\Y,.,} - -P{»,iy,}  (6) 

Bayes 

Plx,\8,.  Y,_,]  - -  plx,\8„  Y.l  O) 

For  output  purposes,  we  can  use  the  lav,  cf  total  probability 

PMY,]^  p[x,\8,=i.  Y,)P{8,  =  H  Y,}.  (8) 


Let  us  take  a  closer  look  at  the  derivation  of  the  above  cycle.  As  u,  and  w, 
are  mutually  independent,  the  Bayes  formula,  which  represents  (6)  and 
(7),  follows  easily  from  (2).  From  the  evolution  of  system  (1)  follows  (5). 
The  Chapman -Kolmogorov  equation  for  the  Markov  chain  8, 

/*{«,-/!  n-.}  =  2  HvP{8,.,-j\Y,-,}  (9) 

J 

which  represents  (3),  can  be  seen  as  a  "mixing.”  To  derive  a 
representation  of  (4)  we  first  introduce  the  following  equation  on  the  basis 
of  the  law  of  total  probability: 

p[*,-,|9,  =  /.  y,-,  1  =  2  [p[Jf,.,|«, .,=>,«,  =  *.  Y,.,  1 

J 

■  P{e,.,=j\8.  =  i.  r,.,}].  no) 

As  8,  is  independent  of  x,  .  ,  if  8,.,  is  known,  we  easily  obtain 
p[x,.,\8,.,=j,  8,  =  i,  Y,.t)=p[x,.l\8,.x=j,  Y,.i J. 

Substitution  of  this  and  of  the  following: 

P{8,.,=j\8,  =  i.  Y,  l)  =  HilP{8,.x=j\Y,.x)/P{6,  =  i\Y,.l) 
in  (10)  yields  the  desired  representation  of  transition  (4) 

P(x, . ,|9,  =  t,  !',.,]  =  £  Ht,P{8l.l~j\Yl.l) 

J 

'  P(x, =A  Y,.l]/P{8l  =  i\Yl.l).  (II) 

Notice  that  the  mixing  of  the  densities  in  (11)  is  explicitly  related  to  the 
above-mentioned  Markov  properties  of  6,  and  the  conditional  indepen¬ 
dence  of  8,  and  x,_ ,,  given  $,.t.  According  to  the  above  filtering  cycle 
there  are  at  any  moment  in  time  N  densities  on  R’  and  /V  scalars.  The 
densities  on  /?"  are  rarely  Gaussian.  Even  if  p(Xo]  Y0]  is  Gaussian,  then 
p(x,|0,  =  /',  Y,]  is  in  general  a  sum  of  A"'1  weighted  Gaussians 
(Gaussian  mixture).  Explicit  recursions  for  these  N'  individual  Gaussians 
and  their  weights  can  simply  be  obtained  from  the  above  filter  cycle. 
Obviously,  the  N  times  increase  of  the  number  of  Gaussians  during  each 
filter  cycle  is  caused  by  (4)  only. 

In  the  sequence  of  ciememai  j  transitions,  (3)  through  (7),  we  can  apply 
a  hypotheses  reduction  either  after  (4),  after  (5),  or  after  (7).  We  review 
these  reduction  timing  possibilities  for  the  fixed  depth  merging  hypotheses 
reduction.  This  fixed  depth  merging  approach  implies  that  the  Gaussian 
hypotheses,  for  which  the  Markov  chain  paths  are  equivalent  durinv  'he 
recent  p-s*  cf  some  fixed  depth,  are  merged  to  one  moment-matched 
Gaussian  hypothesis.  The  degrees  of  freedom  in  applying  this  fixed  depth 
merging  approach  are  the  choice  of  the  depth,  d(  a  1),  and  the  moment  of 
application.  If  the  application  is  immediately  after  each  measurement 
update  pass  (7).  it  yields  the  GPB  {d  +  1)  algorithms  {14],  {16]  In  the 
next  section  we  derive  the  IMM  algorithm  by  applying  the  fixed  depth 
merging  approach  with  depth,  d  =  1 ,  after  each  pass  of  (4).  It  can  easily 
be  verified  that  all  other  timing  possibilities  yield  disguised  versions  of 
IMM  and  GPB  algorithms.  Merging  after  (3)  with  d  =  1  yields  a 
disguised  but  more  complex  IMM  algorithm.  Merging  either  after  (4)  or 
after  (5)  with  d  &  2  yields  a  disguised  but  more  complex  GPBd 
algorithm 


IB.  THE  IMM  ALGORITHM 

The  IMM  algorithm  cycle  consists  of  the  following  four  steps,  of  which 
the  first  three  steps  are  illustrated  in  Fig.  I . 

1)  Starting  with  the  A/  weights  p,(t  -  I),  the  N  means  .?,(/  -  1)  and 
the  N  associated  covariances  /?,(/  -  1),  one  computes  the  mixed  initial 
condition  for  the  filter  matched  to  6,  =  i,  according  to  the  following 
equations: 

A(0  =  5J  >)•  if  A(/)  =  0  prune  hypothesis  i,  (12) 

) 

!)  =  £  l)f;((- D/p.O).  (13) 

l 

)?•(/-  I)=£  //;,£,(/-  !)(*//-  l)  +  (i,(f-l)-i'(t-l))I.  .)T)/fi,V). 

) 

(14) 

2)  Each  of  the  N  pairs  Sd(t  -  1),  R‘(t  -  1)  is  used  as  input  to  a 
Kalman  filter  matched  to  8,  =  /.  Time-extrapolation  yields,  x,(/),  R,(t), 
and  then,  measurement  updating  yields,  R,(t). 

3)  The  N  weights  p,{t)  ate  updated  from  the  innovations  of  the  N 
Kalman  fillers, 

P.U)  =  c-M0  ■  |ft(/)H  ,/Jexp  {-  1/2  J,r(t)Q, ’<r)J,(t)}  (15) 

with  c  denoting  a  normalizing  constant 

d,(t)=y,--A(/)Jf,(f)  (16) 

Q,(t)  =  h(i)R,U)h  r0)  +  g{i)g  r«).  (17) 

4)  For  output  purpose  only,  x,  and  R,  are  computed  according  to 

08) 

*.  =  £  A(')(£<0  + [*.(')- *,) I  ]r]  ’  (19) 

Only  step  I)  is  typical  for  the  IMM  algorithm.  Specifically,  the  mixing 
represented  by  (13)  and  (14)  and  by  the  interaction  box  in  Fig  1.  cannot 
be  found  in  the  GPB  algorithms.  This  is  the  key  of  the  novel  approach  to 
the  timing  of  fixed  depth  hypotheses  merging  that  yields  the  IMM 
algorithm.  We  give  a  derivation  of  the  key  step  1). 

Application  of  fixed  depth  merging  with  d  =  1  implies  that 

^{x,.,|0,-,  =  /.  Y,.x\~N(i.(t-  1).  R,(t-  I)]. 

Substitution  of  this  in  (1 1)  immediately  vields  (13)  and  (14).  with 
f'(t-l)  ft  £{x,.,|*,  =  /.  Y,.,) 
and 

R‘(t - 1) 

the  associated  covariance.  Finally,  we  introduce  the  approximation, 
£(-*,.,1*,  =  /.  Y,.,]~H{S‘(t- I),  R'(t~  I)) 

which  guarantees  that  all  subsequent  IMM  steps  fit  correctly. 

Remark:  The  IMM  ran  be  approximated  by  the  GPB  I  algorithm  by 
replacing  i,(t  -  I)and/?,(f-  I)  in  step  1)  by  S,.  t  and  R,  Together 
with  (12)  this  approximates  (13)  and  (14)  in  step  1)  by,  x'(t  -  I)  =  T,_ , 
and  R'(t  -  1)  at  R'_t.  These  equations  are  equivalent  to  (13)  and  (14)  if 
each  component  of  H  equals  UN,  which  implies  that  8,  is  a  sequence  of 
mutually  independent  stochastic  variables.  The  latter  is  hardly  ever  the 
case  and  we  conclude  that  the  reduction  of  the  IMM  to  GPB1  leads  to  a 
significant  performance  degradation.  Obviously,  the  computational  loads 
of  IMM  and  GPB'  are  almost  equivalent. 
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r  .g  The  IMM  algorithm 
IV  P cKFORMANCE  OF  THE  IMM  ALGORITHM 

Present!,  a  comparison  of  (he  different  filtering  algorithms  for  systems 
with  Markovian  coefficients  with  respect  to  their  performance  is 
ha'  -pered  by  the  analytical  complexity  of  the  problem  [16],  [15].  Because 
f  'his,  such  comparisons  necessarily  rely  on  Monte  Carlo  simulations  for 
specific  examples.  For  our  simulated  examples  we  used  the  set  of  19  cases 
that  have  been  developed  by  Westwood  [18],  To  make  the  comparison 
more  precise,  we  specify  these  cases  and  summarize  the  observed 
performance  results  In  all  19  cases  both  x ,  and  y,  are  scalar  processes, 
which  satisfy  Jr,  -=  a(8,) jr, _ ,  +  b(8,)w,  +  u(r)  and  y,  -  h(d,)x,  + 
g(8,)u,,  with  8,: fl  =»  {0,  1},  u(l)  =  10.  cos  {2x1/100},  x<,  a  Gaussian 
variable  with  expectation  10  and  variance  10,  P{80  -  1}  =  P{80  =  0} 
=  1/2,  while  Hoo  =  (I  -  l/r0)and//n  =  (1  -  1/r,).  The  parameters 
a,  b,  h,  g  and  the  average  sojourn  times  r0  and  r,  of  these  19  cases  are 
given  in  Table  1 

The  results  of  Westwood  ( 18]  show  that,  in  all  19  cases  the  differences 
in  performance  of  the  GPB2  and  the  GPB3  algorithms  are  negligible, 
while  in  only  seven  cases  (5,  6,  8,  16,  17,  18,  19)  the  differences  in 
performance  of  the  GPB 1  and  the  GPB2  algorithms  ate  negligible.  To  our 
present  comparison  the  other  12  cases  (1.2. 3, 4, 7. 9,  10,  II.  12,  13,  14, 
15)  are  interesting.  For  each  of  these  12  cases  we  simulated  the  GPB1,  the 
GPB2,  and  the  IMM  algorithms  and  ran  Monte  Carlo  simulations, 
consisting  of  100  runs  from  I  =  0  to  I  =  100.  For  simplicity  of 
interpretation  of  the  results  we  used  one  fixed  path  of  8  during  all  runs:  8 
=  0  on  the  time  interval  (0,  30] ,  6  =  I  on  the  interval  (3 1 . 60] ,  and  8-0 
on  the  interval  (61,  100] 

The  results  of  our  simulations  for  the  12  interesting  cases  are  as 
follows.  In  six  cases  (I,  2,  7,  12,  14,  15)  both  the  IMM  and  the  GPB2 
performed  slightly  better  than  the  GPB1,  while  the  IMM  and  the  GPB2 
performed  equally  well.  For  typical  results,  see  Fig.  2.  In  the  other  six 
cases  both  the  IMM  and  the  GPB2  performed  significantly  better  than  the 
GPB1.  For  typical  results  see  Figs.  3  and  4.  Of  these  six  cases  the  IMM 
and  the  GPB2  performed  four  times  equally  well  (cases  3,4,  1 1,  and  13) 
and  two  times  significantly  different  (cases  9  and  10). 

On  the  basis  of  these  simulations  we  can  conclude  that  the  IMM 
performs  almost  as  well  as  the  GPB2,  while  its  computational  load  is 
about  that  of  GPB1 .  We  can  further  differentiate  this  overall  conclusion. 

•  Increasing  the  parameters  r0  and  r,  increases  the  difference  in 
performance  between  GPB1  and  GPB2,  but  not  between  IMM  and  GPB2. 

•  If  a  is  being  switched,  then  the  IMM  performs  as  well  as  the  GPB2, 
while  the  GPB1  sometimes  stays  significantly  behind. 

•  If  the  white  noise  gains,  b  or  g,  are  being  switched,  then  the  IMM 
performs  as  well  as  the  GPB2,  while  the  GPB1  sometimes  stays 
significantly  behind. 

•  If  only  h  is  being  switched,  then  in  some  cases  the  IMM,  and  even 
more  often,  the  GPB1  tend  to  diverge  while  the  GPB2  works  well. 

Another  interesting  question  is  how  the  IMM  compares  to  the  modified 
MM  algorithm  and  the  MGEK  filter.  Apart  from  the  GPB  algorithms, 
Westwood  [18]  also  evaluated  four  more  filters,  the  MM,  the  modified 
MM.  the  MGEK,  and  a  MGEK  with  a  "postprocessor  "  For  the  19  cases 
there  was  only  one  algorithm  that  outperformed  the  GPB1  algorithm  in 
some  cases.  It  was  the  MGEK  filter  in  the  cases  1,  3.  and  4.  He  also  found 
that  the  MGEK  filter  performed  in  these  cases  marginally  or  significantly 
less  good  than  the  GPB2  algorithm  As  the  above  experinrems  showed  that 
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TABLE  I 

The  parameters  of  the  19  Cases  of  Westwood  [18] 


CASE 

h-values 

0-DEPENDENT  VALUES 

# 

To 

7  1 

•(0)  .  ■(!) 

b(0l  M 11 

M0)  Ml) 

9<0)  .9(1) 

1 

40 

20 

.995,990 

10 

1.0 

10 

2 

40 

20 

.995.  990 

10 

■!> 

3 

40 

20 

.995,990 

.1 

1.0 

5.0 

4 

200 

too 

.995,  990 

.1 

10 

5.0 

5 

40 

20 

995  990 

80 

1.0 

1.0 

6 

40 

20 

995.990 

1.0 

1.0 

3 

7 

40 

20 

995.900 

S 

1.0 

2.0 

8 

40 

20 

995.750 

10 

1  0 

6 

9 

40 

20 

995 

20 

1  0.95 

5 

10 

40 

20 

995 

10 

10.  80 

2 

n 

40 

20 

995 

5 

1  0..80 

8 

12 

4 

2 

995 

5 

10.. 80 

8 

13 

200 

100 

995 

5 

1  0.  80 

8 

14 

40 

20 

995 

1.5  0 

1  0 

1  0 

15 

40 

20 

995 

1.0 

1.0 

1.50 

16 

10 

2 

95 

.5 

1  0,0  0 

1  0.2  0 

17 

200 

5 

950,0.0 

1.0 

10 

1.0 

18 

60 

5 

950.12 

1.0 

1.0 

1.0 

19 

10 

2 

.95 

5 

VO 

1  0,400 

Fig.  2.  rms  error  for  case  7,  illustrative  of  the  six  cases  ( 1 .  2.  7.  12,  14,  15)  where  both 
IMM  and  GPB2  perform  slightly  better  than  GPB1 


Fig.  3.  rtns  error  foe  case  3,  illustrative  of  the  four  cases  <3.  4.  II.  13)  where  both  IMM 
and  GPB2  perform  better  than  GPB1,  while  IMM  and  GPB2  perform  equally  well 


for  cases  I,  3,  and  4  the  GPB2  and  the  IMM  algorithm  performed  equally 
well,  one  can  conclude  that  the  MM,  the  modified  MM.  the  MGEK,  the 
MGEK  with  "postprocessor,"  and  the  GPB1  are  in  all  19  cases 
outperformed  by  the  IMM  algorithm. 

On  the  basis  of  these  comparisons  one  can  conclude  that  for  practical 
filtering  applications  with  N  —  2.  the  IMM  algorithm  is  the  best  first 
choice.  As  the  IMM  algorithm  has  been  developed  on  the  basis  of  some 
general  hypotheses  reduction  principles,  which  are  Af-invanam,  one  can 
reasonably  expect  that  this  is  also  true  for  larger  N.  But  it  is  unlikely  that 
the  IMM  performs  in  all  applications  almost  as  good  as  the  exact  filter 
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Fig  4.  mu  error  for  case  9.  illustrative  of  the  two  cases  (9  and  10)  where  IMM 
perforins  better  than  GPBl.  but  slightly  worse  than  GPB2  (in  these  two  cases  only  h 
jumps). 


Therefore,  if  the  IMM  performs  not  well  enough  in  a  particular 
application  one  should  consider  using  a  suitable  GPB  (^2)  or  DE 
algorithm  [14],  or  one  might  try  to  design  a  better  algorithm  by  using 
adaptive  merging  techniques  [16].  The  DE  algorithm  might  possibly  be 
improved  by  the  novel  timing  of  hypotheses  reduction  [1].  If  for  a 
particular  application  the  performance  of  the  selected  algorithm  has  a  too 
high  computational  load,  then  it  is  best  to  try  to  exploit  some  geometrical 
structure  of  the  problem  considered  [2],  [11], 

In  situations  where  estimation  has  to  be  done  outside  some  time-critical 
control  loop,  it  is  usually  preferable  to  use  a  smoothing  algorithm  instead 
of  a  filtering  algorithm  [8],  [14],  [2Ij.  In  view  of  the  above  filtering 
results,  this  sugg-  ts  that  the  ideas  that  underly  the  IMM  algorithm  can  be 
exploited  to  develop  better  smoothing  algorithms. 
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Abstract 

>.n  area  of  current  Interest  is  the  estimation  of 
the  state  of  discrete-time  stochastic  systems  with 
parameters  which  may  switch  among  a  finite  set  of 
values  the  parameter  switching  process  of  interest 
Is  modeled  by  a  class  of  semi-Harkov  chains  This 
class  of  processes  is  useful  in  that  >1  pertains  to 
many  areas  of  interests  such  as  the  failure  detection 
problem,  the  target  tracking  problem,  socio-economic 
problems  and  In  the  problem  of  approximating 
nonlinear  systems  by  a  set  of  linearized  models  It 
is  shown  in  this  paper  how  the  transition 
probabilities,  which  #uvern  the  model  switching  at 
each  time  step,  can  be  inferred  via  the  eva.jation  of 
the  conditional  distribution  of  the  sojourn  time 
Following  this,  a  recursive  state  estimation 
algorithm  for  dynamic  systems  with  noisy  observations 
and  changing  structures,  which  uses  the  conditional 
sojourn  time  distribution,  is  derived 

1.  Introduction 

In  this  paper  we  are  concerned  with  failure 
detection  via  recursive  estimation  of  parameters  in 
discrete-time  dynamic  systems.  The  topic  of  interest 
is  stochastic  systems  with  abruptly  changing 
parameters  he.,  model  jumps.  The  recursive  slate 
estimation  algorithm  for  this  problem  developed  in 
this  paper  provides  the  conditional  model 
probabilities  used  for  detecting  the  change  in  system 
parameters  which  signify  component  failures. 

The  abruptly  changing  parameters,  which  switch 
among  a  finite  set  of  values,  are  modeled  as  a  Harkov 
or  a  semi-Markov  chain  with  known  transition 
statistics  (M2.M3.M5-M8.Gl).  Although  the  idea  of 
semi-Markov  chains  is  appropriate  for  the  model 
concerned,  the  analysis  presented  in  the  above  is 
actually  only  for  Markov  chains  (since  the  transition 
probabilities  were  assumed  fixed  and  the  transitions 
depended  only  on  the  latest  state  -  see  Eq.  (8)  in 
( M2 1 ) .  The  process  considered  in  this  paper  is  of 
the  semi-Markov  type  and  pertains  to  many  areas  of 
Interest.  A  failure  in  a  component  of  a  dynamical 
system  can  be  represented  by  a  sudden  change  in  the 
systems  parameters  (Bj.Sl.Wlj.  Also,  a  repair  to  a 
system  represents  a  change  in  the  parameters  1BS|. 
Other  areas  that  this  class  of  processes  pertains  to 
aru  the  target  tracking  problem  (811.  socio-economic 
problems  (G21  and  the  technique  of  approximating 
grossly  nonlinear  systems  by  a  set  of  linearized 
models  (M1.V1.V2I. 

The  first  treatment  of  estimation  in  a  switching 
environment  was  in  (All  where  the  means  and 
covariances  of  the  process  and  measurement  noises 
experienced  Jumps.  As  Indicated  In  (Cl',,  the  optimum 
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state  estimation  in  a  multiple  model  environment  is  a 
function  of  the  elemental  ("model-matched")  state 
estimates  obtained  via  estimators  tuned  to  all 
possible  parameter  histories  .'hus.  with  time,  the 
estimator  must  keep  track  o(  an  exponentially  growing 
number  of  parameter  history  hypotheses  Even  in  the 
case  of  Markov  switching  the  estimation  algorithm 
requires  exponentially  growing  memory  ( T I .  T2) 
Suboptimal  algorithms  like  the  Cencr  alized 
Pseudo-Bayeslan  Algorithm  iGl'81  1*1.  Cl.  T2)  and  the 
Interacting  Multiple-Model  Algorithm  ( 1  MM )  |G2.  83. 

Br)  are  viable  approaches  lo  obtain  a  real-time 
implementable  estimation  algorithm  These  algorithms 
rely  on  different  byopolhe  merging  techniques  lo 
limit  the  memory  and  computational  requirements  |B1) 

In  IS2.C2I  a  semi-Markov  switching  problem  was 
considered,  but  tne  jumps  were  assumed  lo  be 
perfectly  observed  In  (M9|  an  estimation  scheme  for 
semi-Markov  processes  was  developed  based  upon  the 
dtteclion-eslimatioh  algorithm  (OEA)  This  approach 
is  obtained  by  retaining  a  ce;  lain  number  of  most 
likely  parameter  history  hypotheses  The  estimation 
schemes  based  upon  the  DEA  (which  discards  all  but  a 
number  or  most  likely  history  hypotheses)  and  the  GPB 
or  1MM  (which  u  ;e  hypothesis  merging)  algorithms 
represent  different  philosophies  of  algorithm 
design  We  present  an  example  comparing  the  two 
methods  For  a  particular  state  estimation  problem 
later  in  this  paper. 

The  problem  is  formulated  in  Section  2  In 
Section  3  the  sojourn  time  conditional  trobability 
mass  functions  and  the  conditional  transition 
probabilities  which  we  derived  in  (Ml a),  are  given 
here  for  clarity  and  ease  of  reference.  The 
inclusion  of  Section  1.  the  stale  estimation 
algorithm  which  was  developed  in  (Mlbl.  is  for  tht. 
sake  of  completeness.  In  Section  S  simulations  a^e 
presented.  Preliminary  results  on  this  problem  were 
presented  in  (Mia,  Mlbl 

2.  Formulation  of  the  problem 

The  system  is  modeled  by  the  equations 

x(kl  -  FIM(kl)  x(k-l)  ♦  v ( k - 1 ,  M(k))  2.1) 

z(k)  -  H(M(k)l  x(k |  ♦  w(k.Mtk))  (2.2) 

where  M(k)  denotes  the  model  "at  time  k"  -  in  effect 
during  the  sampling  period  ending  at  k  The  process 
and  measurement  noise  sequences.  v(k)  and  w(k).  are 
white  and  mutual!,'  uncorrclated. 

The  model  at  lime  k  is  assumed  lo  be  among  ihe 
possible  r  models 

Mtk)  e  (1.  .r)  (2  31 

For  example 

F(Mlk).j)  ■  F I  (2.1) 

vlk-l.M(k)-j)  ~  N(uy  Q.)  (2  5) 

i.e  ,  the  structure  of  the  system  and/or  the 
statistics  of  the  noises  might  be  different  from 
model  to  model.  The  mean  of  the  noise  ran 

model  a  maneuver  or  a  system  failure  as  a 
deterministic  input 

The  model  switching  process  to  be  considered  here 
is  of  the  semi-Markov  type  The  process  Is  specified 
by  a  Family  of  transition  matrices  P( .  (  t,  ) 
i  e  .  it  is  a  sojourn-time-dependent  Markov"  ISTOM) 


chain,  which  belongs  lo  the  semi-Harkov  class  The 
specification  of  the  STOH  chain  is  more  closely 
related  to  physical  models  because  it  does  not  have 
the  artificial  restart  of  the  sojourn  time  counting 
of  the  semi  Markov  process  for  virtual  transitions' 
and  can  capture  Important  features  in  many  realistic 
situations 

Tor  the  class  of  semi-Harkov  chains  governing  the 
evolution  of  the  system's  model  considered  here,  we 
need  the  pdf  of  the  sojourn  time  conditioned  on  the 
observations,  to  infer  the  transition  probabilities 
The  conditional  transition  probabilities  based  on 
noisy  observations  of  the  system's  slate  are  obtained 
in  the  next  section 

'A  semi-Markov  (SM)  chain  |Ht.  H2.  fill  is 
charadeii  zed  by  a  fixed  matrix  of  transition 
probabilities  (p  1  and  a  matrix  of  sojourn 

•  j 

lime  probability  density  functions 
l  f,  (  t,1  1;  which  are  functions  of  the 
current  stale  i  as  well  as  the  destination  state  | 
of  the  transition  In  a  SM  chain  first  the 
destination  of  the  jump  is  chosen  according  to 

Ip  1  and  then  the  time  after  which  the  jump 
u 

takes  place  (i.e..  the  sojourn  time)  is  chosen 
according  lo  (f  (t,)!.  In  this  model  the 
process  can  undergo  a  virtual  transition  (i  e  .  jump 
in  place'  if  j»i|;  however,  in  this  case,  the 
sojourn  time  counting  is  still  restarted  even  though 
the  system  has  been  in  state  i  for  some  time 

3.  Sojourn  Time  Probability  Mass  Functions  and 
Conditional  Transition  Probabilities 

The  process  M(k).  k*0.1.  .  which  represents 

the  system  model,  can  exist  in  one  of  r 
possible  states  The  current  probabilities  of 
transition  for  the  STDM  process  (chain)  are  functions 
of  the  sojourn  time  t  and  are  defined  as 

p  (t)  -  P(M(k)«j|M(k-l)*i,Ti(k-l)«r)  (3U 

‘J 

where  r,(k-l)  is  the  sojourn  lime  in  stale  i  at 
lime  k-1  It  is  assumed  that  at  k«0  the  sojourn 
time  (in  whatever  state  the  system  model  is)  is 
r *  1  thus  the  values  T  can  lake  are  from  1 
to  the  maximum,  which  at  time  k-1  is  then  k 

i  el  z(i)  be  a  noisy  measurement  of  the  stale  of 
the  dynamic  system  whose  model  undergoes  transitions 
according  to  the  above  described  STOM  process  8ased 
on  the  available  information 

Zk  *  (  z  (  x  )  )‘;1  the  probability  of  the 
model  process  being  in  slate  i  .  denoted  as 
p  |k)  .  is  defined  as 

p(k)=P(M(k|  =  ilZk)  i“l _ r  (3  21 

The  conditional  pmf  of  the  sojourn  time  in  state 
M(k)  =  i  based  on  the  available  information  Z  at 
time  k  is 

gk(T)  i  P<T,(kJ-rlM(k|-=i.Zk)  -  Ptr^kl-rlMIkl-i.Z1'1) 

■  P(M(k-l)«i,.  .M(k-T*X)“i.H(k-T)i*ilM(kl«i,Zk'1)  (3.3) 
where  the  perfect  knowledge  of  the  stale  M(kl 
allaws  one  to  go  down  lo  one  index  less  in  the 

conditioning,  i.e..  Zk  ' . 

Following  (3.1|  the  conditional  probability  of 
transition  from  i  to  j  at  time  k-1  given  the 
observations  Zk  1  is.  in  terms  of  (3.3). 


p  (k-u  £  P(M(k  I*  j|M(k-l  |*(,Z  ‘I 
u 

*  £  P(M(k|«j|M(k-ll*i.T,(k-ll*T.Zk'‘) 
r*l 

■  P<r,(k-ll-r|M(k-l)-i.Zk'') 

•  I  P  ( T )  8k'(T) 
ra  '1 

Note  that  the  argument  of  p.  defined  in 


gkl k ♦  1 )  ■ 


“  i-i  p  Ik-mj 

_  J],  aJkTml 


(3  I)  is  the  sojourn  time  while  the  argument  of 

p  defined  above  is  the  current  tune 
>] 

The  conditional  probability  mass  function  (3  3) 
of  the  sojourn  time  r  in  state  i  at  time  k  is 
given  by  the  following  expressions 

it  (k-1 ) 

bJUI  -  1  '  ITTioT  b',k'n  IJb| 

r  p  (k-s)  ~|  i-i  U  (k-mj 

e:<s‘  -  [_ 1  -  rikjr  b.(ksl  J  n  b'(k  m| 

S-2.  ,k  13  7] 

•  p  (k-m) 

gMk.ll  ■  |1  11,1k. ml  13X1 

expressions  (3.t>)-(3  8)  are  proven  by  induction  in 
(Mia  |  The  ".nations  a,  and  b,  used  above  are 
defined  below 

The  probability  that  the  process  will  slay  s 
time  steps  m  the  same  stale  i  as  it  is  at  time 
k-s  is.  conditioned  on  the  information  at  k-s. 
given  by  the  expression 

b  ik.s)  =  P(M(k)  =  i.  ,M(k-S*l)  =  llM(k-s)  =  i.Z‘  '} 

=  T  "n'p.h)  8,k*‘ln!  s*l.  .k  (3«)l 

n:I  jtn 

Conditioned  on  the  available  information 
z*‘s  at  time  k-s.  the  joint  probability  of 
the  process  residing  in  the  same  slate  i  for  the 
next  s  lime  steps  is  denoted  as 
ai(k.s)=P(M(kj  =  i.  ,Mf k -s ♦U-ilZ* 

*  £p(M|k)  =  l.  ,H(k-SM)  =  i)Mlk-s)  =  J.Z‘  1  )P  (  Ml  k -s  )  =  )IZ‘  ') 
j-i 

=  b  (k.S  Ip  (k-s  I  «  £P(M(k|«i.  .M(V-s*l)  =  i|M(k-s)  =  J.Z“') 

y 

■P^lk-S  / 

=  b,(  k.S )  p(k-s) 

♦  I  1  il  P(M(k  1  =  1. ^  .M(k-sM)  =  ilMlk-s)  =  J.T|(k-s)  =  n.i1''} 

/ 1  nil 

•g  k's(n)Jp  (k-sl 

=  b,  (k.s)  plk-s) 

*  iPl'p  (n)p  (Dp  (21  p  (s-ll  8  k  Mn)  p  Ik-  I 

y,  L  no  J’  “  -I  ' 

-  b,  (  k  ,s  I  p,(  k  -  s  I 

•  X  f  I  P  Inin  P  IDs  k  t(n)  j  P  (k-s) 

y,  L  n.l  >'  1=1  “  1  -J  ' 
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4.  The  state  estimation  algorithm 

As  indicated  in  Sec.  1,  the  optimal  estimator  for 
linear  systems  with  Markov  model  jumps  requires  an 
exponentially  increasing  memory.  Among  the 
suboptimal  approaches  discussed,  it  appears  that  the 
IMM  is  the  most  cost-errective  in  implementation 
[81]  In  view  of  this,  the  state  estimation  for  a 
linear  system  with  so journ-time-dependenl  transition 
probabilities  is  developed  in  the  sequel  based  on  the 
IMM  approach 

In  this  approach,  at  time  k  the  state 
estimation  is  computed  under  each  possible  model 
hypothesis  using  r  fillers  (for  the  r  possible 
models),  with  each  filter  using  a  different 
combination  of  the  previous  model-conditioned 
estimates  Each  model  transition  probability  Is  a 
known  function  of  the  sojourn  time  given  by  (3-D. 
Each  model  has  a  sojourn  time  T,|k)  in  state  i 
which  is,  however,  not  known  The  filler  has  access 
only  to  the  observations  from  which  the  conditional 
pmf  of  the  sojourn  limp  13  8/  .8)  can  be  obtained. 


this  in  turn  is  to  be  used  in  calculation  of  the 
conditional  transition  probabvilities  13  SI 

To  find  the  conditional  pdf  of  the  state  of  the 
dynamic  system  described  by  (2.11-1231  the  total 
probability  theorem  is  used  as  follows 

plx(k)|Zk]  ■  £  p(x(k)|M(k|.|.z(kl.z“  ')  P(M(k|.||Z‘) 

)=' 

-  £plx(k)|M(k|.j.z(k|.Z*  ‘I  p  (k |  M  II 

i=>  1 

re  ,  r  filters  running  in  parallel  The 
model-conditioned  posterior  pdf  of  the  slate,  can  be 
rewritten  as  (with  the  irrelevant  conditioning  on 
Zl  1  In  the  numerator  omitted) 


pi  k( k  )|M| k ) *  j.Z* 


( 1  21 


p|x(k)IM(k)«j,z(k|,Z 
p(z(k)lM(k)«j.x(k)l 

p|z(k)|M(k).j,Z* 
reflecting  one  cycle  of  the  stale  estimation  filler 
matched  to  model  j  starting  with  the  prior,  which 
is  the  last  term  above  The  total  probability 
theorem  is  now  applied  to  this  prior,  yielding 
plx(k||M(k|.j.Z‘"’| 

r 

*  Splx(k)|M(k).j.M(k-l|-i.Zl''|P(H(k-l)-ilM(k|  =  j.Z,  'l 


where 


Iplx(kl|M(k).j.M(k-tl-i.Zll|  p  ,  (k-llk-l) 

i.-i  ,lJ 


M  31 


»|k|  £  P(M( k  1  * j| Zk)  IS  S| 


and 

P ^(k-llk-l!  S  P(M(k-l|.i|M(k|«j.Zl‘ll  (SS) 

Note  that  Eq  (S  3)  represents  a  Gaussian  mixture 
under  the  typical  Gaussian  assumptions  on  the  noise 
terms  in  Eqs.  (2.1)  and  (2.2).  This  mixture  is  then 
approximated  by  a  single  moment-matched  Gaussian- 

Therefore  it  follows  that  the  input  to  the  filter 
matched  to  model  j,  j«l,  ...r.  is  obtained  from  an 
interaction  of  these  r  fillers  This  interaction 
consists  of  the  mixing  of  the  estimates 

x  (k-llk-l)  according  to  the  weightings 

(probabilities!  p.. (k-llk-l)  The 
rlj 

evaluation  of  the  probabilities  h  (I  and  ((.51  in 
the  STOM  situation,  are  the  key  results  needed 
to  obtain  a  recursive  slate  estimation  algorithm  for 
this  type  of  model  switching.  These  probabilities  are 
shown  below  to  follow  from  the  results  in  Section  3 

Pig.  4.1  describes  the  resulting  Interacting 
Multiple  Model  (IMM)  algorithm,  which  consists  of  r 
interacting  filters  operating  in  parallel  The 
mixing  is  done  at  the  input  of  the  filters  with  the 
probabilities,  detailed  later  in  (1.7),  conditioned 
o  n  Zk~'  . 

One  cycle  of  the  algorithm  consists  of  the 
following: 

Starling  with  the  model-conditioned  estimate 
x  (k-llk-l),  with  associated  covariance 
P'lk-ilk-ll.  one  computes  the  mixed  initial 
condition  for  the  filter  matched  to  M(k|*j  according 
to  (1.31  as  follows 

x0l(k-l|k-l)  -  £  x’lk-llk-llp  (k-llk-l  |  (16| 

I- 1 

From  M.5) 

PUj(k-!|k-l)  -  =-P(M(k  )■  j|M(k-l  l*i,Zk  l)P(M(k-l)*i|Zl  ll 

*  =-  P  (k-t)  p  (k-1)  (1  7 1 

C,  u 


'This  is  the  key  step  of  the  IMM  that  yields  an 
algorithm  with  fixed  (and  modest)  computational 
requirements  using  r  filters  It  yields  performance 
comparable  to  the  Generalized  Pseudo  Bayesian 
algorithm  with  r2  filters  (Bll 


where  the  notations  from  (1  1)  and  (3  5)  were  used 
and 

x‘lk-llk-11  iE(x(k-l||M(k-l|.i.Z'1l  (18) 

is  the  model-conditioned  state  estimate  at  time  k-| 

The  expression  of  p  for  the  STOM  case 
■1 

using  terms  involving  sojourn  time  probabilities  is 
the  one  obtained  in  13.51  The  covariance 
corresponding  to  (16)  is 

P0,|k-l|k-I|  •  £p  .  Ik-lIk-lKP'lk-ilk-l) 

i-1  j 

•  (  x'( k-llk-l  |-x0l| k-llk-l ) ) 

-(x'lk-iik-u-xV-iik-nn  (i  9) 

ihe  estimate  (181  and  covariance  (1.9)  are  used 
us  input  to  a  standard  Kalman  filter  matched  to 
M|k)*j  to  yield  the  model-conditioned  estimate 
x’(k|k)  and  its  covariance  P'lklk) 

The  likelihood  functions  corresponding  to  the  r 
filters  are  computed  as 

A  ,1k)  -  plz(k)lM(k).j.z‘  ‘| 

=  PIz(k)|M(k)«  j.x0,(k-l|k-l).P0l(k-l|k-l)|  (HO) 

where  the  past  data  have  beer,  replaced  by  (16)  and 
(1.8)  according  to  the  key  step  of  the  IMM.  The 
model  probabilities  (1.11  are  updated  as  follows 

p  (k)  -  P(M(k)»j|Z*)  *  (k)X  p  (k-1)  p  (k-1)  (1  11) 

1  J  -.j  'J  1 

where  the  conditional  transition  probabilities, 
p  are  as  given  in  (18) 

Eqs  (17)  and  (1.11)  in  combination  with 
p  are  the  key  results  that  make  possible 

the  stale  estimation  for  a  system  with  sojourn-lime- 
dependent  model  transitions 

finally,  for  output  only,  the  latest  slate 
estimate  and  covariance  are  obtained  according  to 
Eqs  (II)  and  (1.3)  as 

x(k|k)  «  £  x'lk|k)  p  (k)  _  (1121 

j=‘  i 

P(klk)  -  £  p  lk)(PJ[k|k) 
i=i  ‘ 

♦  (x'(k|k)  -  x(k|k)HxJ(k|k)  -  x(k|k))  )  |1  13) 

5.  Simulation  Results 


The  algorithm  developed  in  Sec.  1  using  the 
sojourn  lime  pmf  obtained  In  Sec.  3  is  used  to 
estimate  the  state  of  the  system.  In  the  first 
example  the  results  of  this  STOM-based  IMM  estimation 
scheme  are  compared  with  results  obtained  from  an  IMM 
algorithm  based  upon  a  Markov  model  transition 
assumption.  In  the  second  example  the  STDM-based  IMM 
estimation  scheme  is  compared  to  the 
detection-estimation  algorithm  of  IM91.  It  is 
assumed  that  an  STOM  process  described  In  Sec.  2 
governs  the  switching  between  models.  In  the 
following  T  is  the  sampling  period  and  k  Is  an 
integer  representing  the  number  of  sampling  periods 
since  time  zero. 

Example  1 

The  estimation  of  a  controlled  double  integrator 
system  with  process  and  measurement  noises  is 
considered  with  a  gain  failure.  The  two  possible 
modets  are  given  by  the  following  system  equation 

x'( k*l )  -  £  *  ^  x'lk) 

•  £  °i  u(k)  *  £  y2'2  J  v (k )  i-1, 2  (5  1) 

with  measurement  equation 

z(k)  -  (1  0)  x'lk)  ♦  w(k)  (5  2) 

The  models  differ  in  the  control  gain  parameter  b’ 

The  process  and  aeasuremenl  noises  are  mutually 
uncorrelaled  with  zero  mean  and  variances 
given  by 

El  v(k  I  v(j||  -  1-10'*  6,, 


15  3) 


and 


Clwtkl  '*1)11  ■  6k,  <s  ‘,l 

The  control  gain  parameters  were  cliosen  to  be  b'*2 
and  bJ«  l 

The  transition  probabilities  p (1  ( T )  and 
p^lr)  defined  in  13-11  are  shown  in  Tig 
S-l  Note  that  p  (  t  1  foi  ii'j.  are  given 

■I 

by 

p  (t)  “  1  —  p(t)  (S5) 

>j  11 

Thus  we  see  that  Pn ( r  J  is  initially  5  and 
rises  rapidly  to  99  and  then  decreases  towards  I 
which  is  Its  steady  state  value  We  also  see  that 
p  ( r )  has  a  value  close  to  10  for  this  range 
of  t  and  thus  model  state  two  is  essentially  an 
absorbing  state 

Figs  5-2  through  5-4  present  the  results  of 
100  Monte  Carlo  runs  The  true  system  was  initially 
model  I  for  every  run  and  the  model  transitions 
occurred  according  to  the  probabilities  of  Fig  5-1 
For  simplicity,  since  we  are  mainly  interested  in  the 
estimation  of  the  stale,  and  not  in  the  control 
strategy,  we  set  u(kl-3  for  all  k 

The  Markov  based  IMM  used  for  comparison  utilized 
the  a  priori  average  transition  probabilities 
p  (  t  I .  obtained  by  taking  the  expected 

II 

value  of  the  transition  probabilities  shown  in 
Fig.  S-l.  In  other  words,  the  conditional 
probability  p  from  (3. SI  is  replaced  hy  the  a 
priori  (unconditional)  p  given  below  in  (5.71 
The  probability  of  having  a  sojourn  lime  ri 
equal  to  t  is  the  probability  that  model  i  is  in 
effect  for  t-I  steps,  and  then  a  transition  occurs 
at  step  r. 

P(t,=  t)  -  n  p  (t|  J[l  -  pt  (t|]  (5  61 

Thus  we  get 

P  *  2  p  ( r )  Ptr^r)  i  =  l .2  (S  7al 

"  1st  " 

and 

p  *  I  -  pi  (5  7b I 

•)  " 

Figs  S-2  and  5-3  are  plots  of  the  RMS  error  in 
'x,(kl  and  x2(k)  respectively  From  Fig  5-2  we 
can  see  that  the  STOM-based  IMM  estimator  improves 
the  RMS  error  in  Xj(kl  by  as  much  as  20  percent 
From  Fig.  5-3  we  see  that  the  RMS  error  in  x2(k) 
of  the  STOM-based  IMM  estimator  is  as  low  as  one 
third  the  error  of  the  Markov-based  IMM  scheme  Thus 
the  mean-square  error  improved  by  an  order  of 
magnitude. 

Fig.  5-4  is  a  plot  of  the  average  model 
probability  error  This  is  the  error  in  the  filter's 
determination  of  the  correct  system  model. 

Typical  running  limes  for  the  STOM-based  IMM  vs. 
the  Markov-based  IMM  are  in  the  ratio  of  3:1.  The 
length  of  the  lime-span  over  which  the  sojourn  lime 
pmf  is  computed  can  be  truncated  -  it  becomes 
negligible  after  15  steps  This  keeps  within 
reasonable  limits  the  additional  calculations  of  the 
STOM-based  filler  and  prevents  any  growth  of  the 
computational  or  memory  requirements. 

Example  2 

In  this  example  we  make  a  comparison  between  the 
detection-estimation  algorithm,  (OEA),  based 
semi-Markov  estimator  of  |M9)  with  the  STOM-based  IMM 
estimator  of  this  paper.  For  this  purpose  the  system 
and  the  semi-Markov  model  switching  process 
attributes  are  as  in  (M9|  example  3.  and  are  repeated 
here  for  ease  of  reference. 

The  model  process  M(k)  is  taken  as  a  semi-Markov 
chain  The  scalar  system  is  described  by  (M9| 

x (k ♦  I )  -  I  04  x(k)  ♦  v ( k  1 

z(k  |  •  x(k )  ♦  0 ( M ( k |  )w(k  |,  K-0,1.2,  (5  81 

where  r  •  3  models,  0(11-100,  0(21-10.  and  0(31-1 


Here  ( v(k ) }  and  (w(k|)  are  mutually  independent 
zero-mean  Gaussian  white  noise  sequences  with 
covariances  Q-0  1  and  R-l  0,  respectively  The 
initial  conditions  are  x(Q)~f/(30.4G0).  P(M(0)-l)*l/3 
for  i » 1.2.3  For  the  real  system  x(0)*l  in  every 
simulation  The  process  M(k|  is  modeled  by  a 
semi-Markov  chain  with  the  imbedded  Markov  chain 
transition  probabilities  given  byPn*Pjj"P33-0.  P|2'0 .7. 
Pn«0.3,  p^-0.6.  pjj-0.4.p]i*0  3,  and  p^-0  7  The  sojourn 
time  probability  mass  functions  p^t)  are  assumed 
to  be 

P^t)  *  a,expl  - 1  t  -  3 1 1 
P^t)  =  a2exp|-|T-6l| 

Pj( t I  -  a3exp( -It -81 1  (5  9 1 

for  riO  with  a  such  that 

Y  p  (  t  1-1.  !*1.2.3  IS  1 0 1 

■-a  1 

The  results  of  SO  Monte  Carlo  runs  average  are 
shown  in  Tigs  5-5.  5-6  In  Fig  5-5  we  compare  the 
rms  state  errors  of  the  two  filter  OEA  based 
semi-Markov  estimator  of  |M9I  with  our  two  filler  GP8 
based  semi-Markov  approach,  and  with  the  GPB 
estimator  using  3  filters  Note  that  the  values  for 
the  OEA  estimator  are  two-time-step  smoothed  values 
(see  [M9 1.  Fig  7.  M-2  most  likely  histories 
retained)  whereas  the  values  for  the  STQM-IHM 
estimator  are  filtered  values  We  can  see  that  our 
estimator  with  two  filters  is  stable  as  opposed  to 
the  unstable  two-filter  DEA  method 

The  plot  of  the  3  filter  ST0M-IMM  estimator  shown 
in  Fig  5-S  is  given  so  that  one  can  compare  the 
improvement  obtainable  by  adding  an  extra  filter  to 
this  approach.  We  see  that  the  long  term  trend  is 
f or  the  3  filter  STDM-IMM  to  givp  a  smaller  rms  error 
than  the  version  with  2  filters. 

In  Fig  5-6  we  compare  the  probability  of  error 
obtained  using  a  4  filter  OEA  estimator  versus  the  3 
filler  STDM-IMM  estimator.  Both  curves  were  obtained 
from  a  filtering  operation  (see  [M9I  Fig.  10.  N  =  0). 

We  can  see  that  the  present  estimator  gives  a  much 
clearer  indication  of  the  correct  system  structure 
and  hence  is  preferable  for  failure  detection 
G.  Conclusion 

We  have  applied  the  recursive  slate  estimation 
algorithm  for  dynamic  systems,  whose  stale  mcdel 
experiences  jumps  according  to  a  sojourn-lime- 
dependent  Markov,  STDM.  chain,  to  the  problem  of 
failure  detection.  The  algorithm,  which  is  of  the  IMM 
type,  uses  noisy  state  observations  and  the 
calculations  are  done  in  the  following  order: 

1.  Probability  of  each  model  being  the  current 

model 

2.  Sojourn  time  pmf  in  the  current  model 

3.  Model-conditioned  state  vector  estimates  and 

covariances 

4.  Overall  state  vector  estimate  and  Its 

covariance. 

The  first  example  simulated  indicates  that  the 
use  of  the  STOM-based  IMM  estimator  can  give  a 
substantial  improvement  in  slate  estimation  over  a 
Markov-based  IMM.  The  latter  relies  on  the  a  priori 
average  transition  probabilities  while  the  former 
uses  conditional  transition  probabilities  obtained 
from  the  conditional  sojourn  time  distribution.  This 
example  shows  that  the  STOM-based  scheme  Is 
substantially  belter  than  the  Markov-based  scheme  In 
determining  the  true  system  model,  which  Is 
beneficial  for  failure  detection  schemes 

The  second  example  simulated  shows  that,  for  the 
particular  system  under  consideration  the  STOM-based 


IMM  estimator,  which  is  an  hypothesis  merging 
technique,  compares  favorably  in  terms  of  the 
probability  of  error,  to  the  detection-estimation 
algorithm  based  estimator,  which  discards  the 
unlikely  parameter  history  hypothesis 
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Figure  S-S  RMS  slate  errors  *0'  example  2  Note  that 
m«2.  stands  for  a  two-filler  STON-IMM  based  estimator  etc 


0.00  zoo  i  .oo  6.i 


HH  1H.IM  IZ  ert  H.tJfl  u.MUl.wtfl.ea 


S  00  10.0«  IS  00  20.00  2S.00  30 


Figure  s-2  RMS  error  in  x,(k|.  STOM-based  Ihh  and  Narkov-based 
Imm  Markov-based  IMM  errors  are  shown  by  the  solid  line 


Itgure  S-t>  Probability  of  error  m  system  structure  detection 
for  example  2  Note  lhai  UCA4  stands  for  the  ^-filter  OTA 
based  estimator  etc 
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Figure  S-3  RMS  error  in  x,(k|.  STOH-based  IMM  and  Markov-based 
IMM  Markov-based  IMM  errors  are  shown  by  the  solid  line 
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Figure  S-4.  Average  model  probability  error  magnitudes. 
STOM-based  IMM  and  Markov'based  IMM  Markov  based  IHH 
errors  are  shown  by  the  solid  line 
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Distributed  Adaptive  Estimation  with 
Probabilistic  Data  Association* 


K.  C.  CHANGt  and  Y.  BAR-SHALOM$§ 

A  fusion  algorithm  for  target  state  estimation  under  cluttered  environment 
with  uncertain  measurement  origins  and  uncertain  system  models  in  a 
distributed  manner  can  be  applied  for  tracking  a  maneuvering  target  in  a 
cluttered  and  low  detection  environment 


Key  Words — Distributed  estimation;  multiple  model;  target  tracking;  probabilistic  data  association; 
Bayesian  methods;  distributed  sensor  networks. 


Abstract — The  probabilistic  data  association  filter  (PDAF) 
estimates  the  state  of  a  target  in  a  cluttered  environment. 
This  suboptimal  Bayesian  approach  assumes  that  the  exact 
target  and  measurement  models  are  known.  However,  in 
most  practical  applications,  there  are  difficulties  in  obtaining 
an  exact  mathematical  model  of  the  physical  process.  In  this 
paper,  the  problem  of  estimating  target  states  with  uncertain 
measurement  origins  and  uncertain  system  models  in  a 
distributed  manner  is  considered.  First,  a  scheme  is  described 
for  local  processing,  then  the  fusion  algorithm  which 
combines  the  local  processed  results  into  a  global  one  is 
derived.  The  algorithm  can  be  applied  for  tracking  a 
maneuvering  target  in  a  cluttered  and  low  detection 
environment  with  a  distributed  sensor  network. 

1.  INTRODUCTION 

The  major  difficulty  in  tracking  a  target  with 
switching  models/parameters  in  a  cluttered 
environment  is  due  to  the  fundamental  conflict 
between  the  operations  of  model/parameter 
identification  and  data  association,  since  the 
measurements  with  large  innovations  are  con¬ 
sidered  as  unlikely  to  have  originated  from  the 
target  of  interest.  In  this  paper,  a  multiple  model 
approach  in  conjunction  with  the  probabilistic 
data  association  (PDA)  filter  (Bar-Shalom  and 
Tse,  1975;  Bar-Shalom,  1978)  to  track  a  target 
with  switching  models  using  distributed  sensors, 
is  presented. 
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Several  approaches  have  been  proposed  to 
perform  the  state  estimation  of  a  system  together 
with  identification  of  each  model  (out  of  a  finite 
set)  in  a  centralized  framework.  One  of  the 
significant  schemes  is  the  so-called  generalized 
pseudo  Bayes  (GPB)  method  (Tugnait,  1982; 
Chang  and  Athans,  1978)  and  the  other  is  the 
interacting  multiple  model  (IMM)  algorithm 
(Blom,  1984;  Blom  and  Bar-Shalom,  1988).  The 
general  structure  of  these  algorithms  consists  of 
a  bank  of  filters  for  the  state  cooperating  with  a 
filter  for  the  parameters.  A  GPB  algorithm  of 
order  n  (GPBn)  needs  AT  filters  in  its  bank 
(Tugnait,  1982).  The  IMM  algorithm  performs 
nearly  as  well  as  the  GPB2  method  with  notably 
less  computation,  namely,  at  the  cost  of  GPB1 
(Blom  and  Bar-Shalom,  1988).  A  distributed 
estimation  scheme  with  uncertain  models  has 
also  been  derived  (Chang  and  Bar-Shalom, 
1987).  However,  in  all  the  above  approaches,  a 
perfect  data  association  was  assumed,  i.e.  there 
is  no  uncertainty  in  measurement  origins. 

To  take  into  account  the  data  association 
problem,  an  adaptive  PDA  algorithm  was 
presented  in  Gauvrit  (1984)  for  tracking  in  a 
cluttered  environment  with  unknown  noise 
statistics.  This  algorithm  identifies  on  line  the 
unknown  variances  of  the  process  and  measure¬ 
ment  noises  but  uses  an  earlier  (static)  multiple 
model  approach  (Bar-Shalom,  1988).  In  this 
paper,  a  distributed  estimation  problem  which 
takes  into  account  both  model  and  measurement 
origin  uncertainties  will  be  derived.  To  handle 
the  model  uncertainty,  a  more  general  formu¬ 
lation  with  dynamic  multiple  models  described  by 
Markovian  parameters  will  be  adopted.  These 
parameters  may  switch  within  a  finite  set  of 
values  which  represent  different  system  models. 
To  take  care  of  the  missing  and  false 
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measurements,  the  PDA  scheme  will  be 
employed.  The  probabilities  of  associating 
measurements  to  a  target  given  different  system 
models  will  be  computed  and  used  to  weight  the 
combination  of  state  estimates. 

The  problem  is  formulated  in  Section  2.  A 
centralized  algorithm  which  combines  the  IMM 
algorithm  and  the  PDA  filter,  resulting  in  the 
MMPDA  (multiple  model  PDA)  filter,  for  local 
processing  will  be  described  in  Section  3.*  Then 
the  fusion  algorithm  which  combines  the  local 
processed  results  from  multiple  sensors  into  a 
global  one  will  be  presented  in  Section  4. 

The  algorithm  can  be  applied  for  tracking  a 
maneuvering  target  in  a  cluttered  and  low 
detection  environment  with  a  distributed  sensor 
network  (DSN). 

2.  PROBLEM  FORMULATION 

Let  us  consider  the  two-node  scenario  similar 
to  that  given  in  Chang  et  al.  (1986),  where  each 
node  processes  the  local  measurements  from  its 
own  sensor  and  sends  the  local  estimates  to  the 
fusion  processor  periodically.  The  fusion  pro¬ 
cessor  then  sends  back  the  processed  results 
after  each  communication  time. 

The  dynamics  of  the  target  in  track  are 
modeled  as 

x(k)  =f[x(k  -  1),  M(k),  v[M(k),  k  -  1]]  (1) 

where  x (k)  is  the  state  vector,  v[M(k),  k  -  1] 
the  process  noise  vector  and  M{k)  the  system 
model  from  time  k  -  1  to  k.  Assume  the  random 
model  process  M(k)  is  Markov  and  it  can  only 
take  values  from  a  finite  set  M,  which  contains  r 
distinct  models,!  i.e. 

M  -{My};.,.  (2) 

The  measurement  system  is  modeled  as  follows. 
If  the  measurement  originates  from  the  target  in 
track,  then 

z‘(k)  =  h‘[x(k),  M(k)]  +  w'[M{k),  &]  (3) 

where  z‘(k)  is  the  measurement  vector  from 
sensor  i  and  tv^A/^),  fc]  is  the  corresponding 
measurement  noise  vector.  The  two  noise 
sequences  are  mutually  independent  and  inde¬ 
pendent  of  the  initial  state. 


*  The  MMPDA  algorithm  has  been  implemented  in  the 
interactive  software  MULTIDAT  (Bar-Shalom,  1987,  1988). 

tThe  models  can  have  states  of  different  dimension.  In 
this  case,  the  lower  dimension  state  vectors  are  augmented 
with  suitable  components  that  are  zero  w.p.l,  to  make  them 
compatible.  This  is  elaborated  on  in  Section  5. 

t  Such  a  rule,  also  called  “gating",  considers  only  the 
measurements  within  some  distance  from  the  predicted 
measurements  (for  details,  see,  e  g.  Bar-Shalom  and 
Fortmann  (1988)). 


As  in  the  PDA  filter,  it  is  assumed  that  a  rule 
of  validation  of  the  candidate  measurements!  is 
available  such  that  it  guarantees  that  the  current 
return  will  be  retained  with  a  given  probability. 
For  each  sensor,  denote  the  validated  measure¬ 
ments  at  time  k  as 

Z'(k)=  {z;(k)};mj,  (4) 

where  m‘k  is  the  number  of  validated  measure¬ 

ments  of  sensor  i  at  time  k,  and 

Z'*  —  {Z'(/)}*=,-  (5) 

The  local  model-conditioned  state  pdfs  at 

sensor  i  are 

P(x(k)  |  Mj(k),  Z‘  \  Y‘  k), 

i  =  l,2;  ;  =  1 . r  (6) 

with  the  corresponding  model  probabilities 

P{M,(k)\Z‘k,  Y‘k}, 

(  =  1,2;  /  =  1 - r  (7) 

where 

r*  =  {Y'(l) _ >"(*)}  (8) 

and  Y'(k)  denotes  the  information  received  by 

node  i  during  the  sampling  period  ending  at  time 
k,  which  is  defined  as  the  fusion  result  (namely, 
global  conditional  pdf)  up  to  time  k  —  \. 

Assuming  lossless  communication  and  that  the 
information  communicat'd  is  the  sufficient 

statistics,  i.e.  the  information  contained  in  Y‘,k  is 
equivalent  to  the  information  in  Z'  *~\  then  we 
have  the  following  equality: 

p(x(k)  |  Z'  *-\  Y‘  k)  =  p(x(k)  |  Zi  k~\  Zlk~') 
~P(x(k)  |  Z*-1)  (9) 

where  i  represents  all  sensors  other  than  sensor  z 
and  Z*  =  {Z(/)}/*«i,  where  Z(/)  represents 
measurements  from  all  sensors  at  time  /. 

Given  the  above  models,  the  question  now  is 
how  the  global  conditional  pdf  can  be  con¬ 
structed  by  fusing  together  the  local  ones. 
Specifically,  we  shall  investigate  what  is  the 
necessary  and  sufficient  information  that  has  to 
be  transmitted  between  nodes.  The  derivations 
will  be  carried  out  for  arbitrary  pdfs;  however, 
the  simulations  assume  linear  models  with 
Gaussian  random  variables,  in  which  case  the 
state’s  model-conditioned  pdf  (6)  is  Gaussian 
and  the  overall  conditional  pdf  of  the  state  is  a 
Gaussian  mixture  (Bar-Shalom,  1988). 

3.  CENTRALIZED  ALGORITHM  FOR  LOCAL 
PROCESSING 

For  each  local  node,  the  centralized  algor¬ 
ithm  where  all  measurements  are  sent  to  and 
processed  with  one  processor  is  described  below. 
The  goal  is  to  compute  the  conditional  state 


Distributed  adaptive  estimation 


361 


distribution  given  the  local  accumulated  measure¬ 
ments.  With  only  model  uncertainty,  the  local 
conditional  pdf  at  sensor  i  can  be  obtained  as 

p(x(k)  |  Z'*,  r  *) 

=  2  P(x(k)  |  Z‘\  Yi  k) 

i- 1 

x  P{M,(k)\Zik,Yik}.  (10) 


When  the  additional  measurement  origin  uncer¬ 
tainties  are  present,  the  above  equation  becomes 

p(x(k)  |  Zik,  Y‘  k) 

=  £{2p(*(*)|M/(*).  0',,Z'k.  Y,k) 

1= 1  1 01. 

x  P{ei\Mj(k),  Z‘k,  Y'k } J 
x  P{Mi(k)\Z,  k,Y,  k}  (11) 


where  &',  is  the  event  that  z[(k)  is  the  correct 
measurement  and  8‘0  denotes  no  correct 
measurement. 

The  first  term  on  the  right-hand  side  of 
equation  (11)  is  the  standard  PDA  filter  based 
on  model  Mr  where  for  each  8 \ 


p(x(k)  |  Mj(k),  8{,  Z,  k,  Y,  k) 


1 


p[Z\k)\x[k), 


Mj(k),  81  Zi  k~\  Yi  k) 
xP(x(k)  |  Mi(k),  Z''k~\  Yl  k) 


(12) 


where  8\  has  been  omitted  in  the  last  term  above 
(since  it  is  irrelevant)  and 

cW(k),  8',] 

=  fp(Z‘(k)  |  x(k),  Mj(k),  8\,  Zi  k~\  Y‘  k) 

xp(x(k)  |  Mj(k),  Zi  k~\  Yi  k)  dx(k) 

=  p(Z'(k)  |  Mj(k),  8[,  Zik~\  Y'-k).  (13) 


Using  Bayes’  rule,  the  second  term  on  the 
right-hand  side  of  equation  (11)  is 


P{8[  |  Mj(k),  Zl  k,  F*} 
p(Z‘(k)  |  M/ik),  81  Z'  k-\  Yi  k)P{8‘,  |  M,(k), 
= _ Zi,k~\  Y,  k)p(Mi{k),  Z*-*-,  Yik) 

P(Z'{k)  I  M,(k),  Zik~\  r  k) 

xpiMfik),  Z'  *-,  Yi  k) 


1 


p(Z‘(k)  I  M'ik),  8lZi  k~\  Y‘  k) 
!, !  M£k).  Z’  ‘ 

c‘\[Mj(k),  8\\ 


c'2[M,{k)\ 

x  P{8[  j  Mjik),  Z'  k~\  Y,  k) 
1 

c'2(  Af,(/c)) 


x  P(8,li\  M,(k),  Z,k-\  Y'  *}  (14) 


where 

c'2[M,(k)]  =  2  c',[A^(Ac),  8{] 

01, 

xP{8[\Mi(k).  Zik~\  Yik} 

-  P(Z‘(k)  |  Mj(k),  Z'  *_1,  Y,  k).  (15) 


In  equation  (13),  the  joint  measurement  density 
is  (see,  e.g.  Bar-Shalom  (1988)) 


p{Z‘{k)\Ml{k)h  8[,  Z'k-\  Yik) 

~  fl  P(z'i(k)  |  M^k),  8[,  Z'  k~l,  Y'  k) 

/=  1 

=  f  v r'  if  /,  -  0 

1  VI m'i*'p[z'li(k)  |  A/,(/c)]  otherwise 


where  Vk  is  the  volume  of  the  validation  region, 
because  our  assumption  on  the  incorrect 
measurements  being  uniformly  distributed,* 
independent  from  each  other  and  from  the 
correct  measurement,  and 

p[z[(k)  |  Mj(k) ] 

~Pcp{zm\M,{k),  8\,  Z‘k~\  Y‘k)  (17) 


is  the  truncated  density  which  is  zero  outside  the 
validation  region  where  PG  is  the  probability  that 
the  correct  return  will  lie  in  the  validation 
region. 

In  equation  (14),  P{6‘, J  Mf(k),  Z‘  k~l,  Y'  k}  is 
the  prior  probability  of  the  event  based  on 
model  M,  to  be  correct  at  time  k.  By  choosing  a 
large  enough  validation  threshold,  this  prob¬ 
ability  becomes  independent  of  Mt(k)  and  is 
assumed  to  be  the  same  for  all  0},  unless  target 
signature  information  can  be  used.  If  no  such 
information  is  available,  then 


P{di\Mi(k),  Zik'\  Yik) 


if  /,  =  0 
otherwise 


(18) 


where  PD  is  the  probability  that  the  correct 
return  will  be  detected. 

For  each  model  Mt{k)  and  event  8[,  equation 
(12)  is  the  standard  filtering  equation.  In  that 
equation,  by  using  the  IMM  approach  (Blom 
and  Bar-Shalom,  1988),  the  extrapolated  pdf  is 
obtained  by  combining  the  extrapolations  of  the 


*  For  more  elaborate  models  see  Bar-Shalom  (1988). 
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prior  pdfs  (independent  of  the  event  9)) 
p(x(k)  |  Mj(k),  Zi  k~\  Y‘k) 

-  t  P(x{k)  |  M,{k),  M,(k  -  1),  Z'  *-',  Y'  k) 

/-  1 

x  P(M,(k  -  1)  |  M,(k),  Z‘-k-\  Y,k} 

i,  ''(*)  I  M,{k),  M,(k  -  1),  Z‘  k~\  r  k) 

x  P{M,(k),  M,{k  -  1)  |  Z‘k~\  Y‘-k) 
P{M,(k)  I  Z‘  k-\  >"*} 

-ra!w‘wl'w- 

M, (k  -  1),  Zik~\  Yik) 

<  P{M,(k )  |  M,(k  -  1)} 
xP{M,(k-\)\Zik-\Yik}\  (19) 

where  p(x{k)  |  Mt(k),  M,(k  -  1),  Z‘-k~\  Yi  k)  is 
the  extrapolation  of  the  conditional  state  pdf 
given  Z‘  k~l  and  Y‘  k  from  model  M,(k  —  1)  to 
model  Mj{k)  and 

cii[Mi{k))  =  P{Mi{k)\Zik-\  Yik) 

=  X  P{M,{k)\M,{k  -  1)} 

/=1 

y.r{LMl(,k-\)\Z,-k-\  y-*}.  (20) 

The  last  term  of  equation  (11)  is  the  a 
posteriori  model  probability,  which  is  obtained 
as 

p[M,(k)  | zi  k~\  r  *> 

=  ^p(Z‘(k)\Ml(k),Z‘  k-\  r‘) 

c* 

xP{Mj(k)  |  Zi  k~\  Yi  k } 

=  ~ci[M;(A:)]cUM/(A:)]  (21) 

c* 


Fir,.  1.  Centralized  MMPDA  algorithm  with  r  =  2  at 
sensor  i 


where 

=  p(Z'(*)|Z'*-\  r‘)  (22) 

and  CjfA/^/:)]  and  c^A/^/:)]  have  been  obtained 
in  equations  (15)  and  (20),  respectively. 
Equations  (12)— (21)  complete  a  recursive  cycle 
of  the  local  processing.  A  flow  diagram  of  the 
local  MMPDA  algorithm  is  given  in  Fig.  1.  The 
flow  of  data  is  represented  by  the  model- 
conditioned  means  i,  and  the  model  prob¬ 
abilities  Pr 

4.  FUSION  ALGORITHM 

With  the  local  conditional  pdfs  obtained  in 
Section  3,  we  can  now  derive  the  fusion 
algorithm  to  obtain  global  pdf.  Similar  to 
equations  (10)  and  (11),  the  global  conditional 
pdf  can  be  obtained  as 

p{x(k)  |  Zk) 

=  2  p(x(k)\  M,(k),  Zk)P{M,(k)\  Zk) 

/=! 

=  2\22p(x(k)\M,(k),  el  el  zk) 

/-i  le),  efj 

xP{elel\Mi(k),Zk}}p{Mi(k)\Zk} 

(23) 

Assuming  measurements  from  different  sensors 
are  independent  given  the  target  state,  then  the 
first  term  on  the  right-hand  side  of  equation  (23) 
can  be  obtained  as 

p(x(k)  I  el  el  Zk) 

_ l 

~c[Mi(k),  el  el) 

xp(Z(k)\x(k),M,{k),  el  el  zk~l) 
xP(X(k)\Mi(k),  el  el  z*-) 

_ i 

~c[M,(k),  el  el) 

xfUp&WlxW,  Miik),  61  Z*-1)] 

l-l 

xp{x{k)  |  Z*_1) 

1 

“  cWk),  el  el) 
n[/7(zw|x(fc),Af,(fc),  eizk~l) 

_ xP{x(k)\M,{k),  Z*-1)} 

x  p{x(k)  |  Mj(k),  Z*"1) 

(24) 
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where 

c[Mj(k),  el  el\ 

=  J p(Z(k)\x(k),  M,(k),  el  el  Zk~') 

Xp(x(k)  |  Mjik),  el  el  Zk~l)  dx(k) 

=  p(Z(*)|  *#,-(*),  eJ.X.Z*"')  (25) 

is  the  normalization  constant. 

Since  from  equations  (12)  and  (9) 

p(x(k)  I  M,{k),  e\,  Z'-‘,  Yik) 

1 

~c\[Mi{k),  e;j 

xP(Z‘(k)\x(k),Mi(k),  0i  Zk~‘) 
Xp(x(k)m(k),Zk-‘).  (26) 

Equation  (24)  can  be  rewritten  as 
P(x(k)\M,(k),  el  e\,zk) 
i 

_c[a/;(*),  el  ej2] 

n  (c^ik),  0',] 


xp(x(k)\M,(k),  e[,zik,  Y‘  k)\ 

p(x(k)  I  Mj(k),  Zk~') 


1 


(27) 


co[Mj(k),  el  el 
n  P(x(k)  |  M,(k),  e{,  z'  \  yi  k) 

i*  l 

P(x(k)  |  M,(k),  Zk~l) 
where  the  denominator  can  be  derived  as 

p(x(k)  |  M,(k),  Zk-‘)  =P(*(k)>MjW\  zk  *) 
1/1  '  p(Mj(k)  |  Z*-1) 

t  p(x(k)  j  Mj(k),  M,(k  -  1),  Z*'1) 

i-1 

_  x  P{M,(k )  j  -  l)}P{Af,(fc  -  1)  1  Z*"1} 
i  P{M,(k)  I  M,(k  -  1)}P{M,(A:  -  1) !  Zk~1} 

(28) 


and 

CoWi(k),  el  0,2J  = 


,  n;,  _  c(Mt(k),  el  el 


n  clM^k),  6‘] 


=1- 


n  p(x(k)  I  M,(k),  e\,  Z'  k,  Y(k) 


p(x(k)  |  M,(k),  Zk~’) 


d x(k) 


(29) 


Assuming  0't>  and  8f2  are  independent  given 
the  target  state,  then  similarly  to  Chang  et  al. 
(1986),  the  second  term  of  equation  (23)  can  be 
obtained  as 

P{dlel  \M,(k),zk) 

~7i^m  lpW-  6‘-  Z(*>  |x(*)' mi 

Zk~‘)p(x(k)  |  Af((&),  Zk~')dx(k) 

\ 

~c,[K(k)l 

hP(x(k),eiZ‘(k)\M,(k),  zk-') 

XJ  p(x(k)\M,(k),Zk-')  <k(*) 

np(x(k)\M;(k),  e\,  Z’(k), zk-') 

xi  dJf(A) 

(30) 

where 


ci(H(^^  ~ P(z(k)  |  M,(k),  Zk~l)  (31) 

and 


c2[a/,(*)]=- - - : 

n  p(z'(k)\M,(k).zk-‘) 

/= i 

n  cifw^/c)] 

i  =  i 


(32) 


are  normalization  constants,  where  Cj[Mi(k)] 
was  given  in  equation  (15). 

Since  the  information  contained  in  Z*~*  is  the 
same  as  that  in  {Zik~\  Yik)  (see  equation  (9) 
for  details),  equation  (30)  can  be  written  as 


P{8l  81 1  Mj(k),  Zk) 

n  P(e[\Mi(k),  z,k,  r-*} 

i- 1 

=  c2[M,(k)\ 

.  U  p(x(k)  |  M,(k),  0j,  Z‘  k,  Y‘  k) 

XJ  p(x(k)\M,(k),  Zk  l)  dr(*) 

x  c„[m, (*).«;,.  o;j  (33) 


is  the  new  normalization  constant. 


From  equations  (27)  and  (33),  equation  (23)  can 
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be  written  as 


r 

p(x{k)  |Z*)  =  £ 
1 


1 

c2[M,(k)\ 


n  [p(x{k)\M,(k),  0}..Z*-\  Y,k)  ' 

1  =  1 

x  P{d[\M,(k),  Z,k,  Y‘k) I  1 
p(x(k)  \M,(k),  2‘-') 


xP{M,(k)\Zk). 


(34) 


The  last  term  of  equation  (34)  is  the  global  a 
posteriori  model  probabilities.  With  equations 
(31)  and  (32)  we  have 

P{M,(k)\Zk) 

~zp{Z.{k)  |  M,(k),  Zk~')P  {Mfk)  |  Zk~') 
c 


=  '-ct[M,(k))P{Ml(k)\Zk-') 


c2\M,{k)} 

c 

I  *=  l 

'  P(M,(k)  JZ1^} 

c2M(*)1 

c 


n  [p(Z‘(k)\  Mf(k),  Z'  k~\  Y'  k) 

xP{M,(k)\Zk-'}\ 

P{M,{k)\Zk-') 

c 

n  [  p{M,{k)\zi(k),  zk-'} 

_ xp(Z'(fc)|Z‘-‘)j 

P{M,(k)  |  Z*'1} 

c2{M,(k)\ 

c 

n  [P{Mj(k)  |  Z‘  k,  Y,  k}p(Z‘(k)  |  Z*-‘)J 
1-1 

P{M,(k)\Zk~'} 

{U(,  u  I!  P{M,(k)\Z‘k,  Y,k) 

cz(M,  (*)]-» _ n, 

c'  P{M,(k)\Zk~'} 


where  the  denominator  is  the  same  as  that  of 
equation  (28)  and  the  normalization  constants  c 


and  c’  are 

c  =  p(Z(k)  |  Z*~*) 
and 

c  c 

c'  =  — - =  — - 

n  p(zw|z‘-)  nc; 

i-i  i  - 1 

4. 1 .  Overview  of  the  fusion  algorithm 

From  the  above,  it  follows  that  the  global  a 
posteriori  pdf  and  model  probabilities  are 
obtained  by  combining  (multiplying)  the  local  a 
posteriori  pdfs  and  model  probabilities  and 
removing  (dividing)  the  common  a  priori  pdf  and 
model  probabilities.  From  equation  (34),  we  can 
see  that  for  each  model,  the  conditional  global 
pdf  given  that  this  model  is  correct  is  obtained 
by  the  sum  of  global  fused  pdfs  given  all  possible 
global  event  pairs  6}lt  Q\.  The  overall  global  a 
posteriori  pdf  is  then  obtained  by  the  sum  of 
global  pdfs  of  each  model  weighted  by  the  global 
a  posteriori  model  probabilities.  Equations  (34) 
and  (35)  represent  the  complete  cycle  of  fusion 
processing.  From  them  it  follows  that  the 
information  needed  to  be  communicated  from 
local  nodes  to  the  fusion  node  consists  of: 

(a)  the  model  probabilities; 

(b)  the  association  event  probabilities;  and 

(c)  the  corresponding  pdfs  (mean  and  covari¬ 
ance  for  Gaussian  case). 

A  summary  flow  diagram  of  the  fusion 
algorithm  with  two  models  is  given  in  Fig.  2.  For 


(36) 

(37) 


«  Ih  1|k  II  [k  U  k  1  I  P  Ik  tj  h  \  1 


Fig.  2  Distributed  MMPDA  algorithm  with  r  -  2 
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simplicity,  only  the  mean  of  each  pdf  is  shown  in 
the  figure.  References  to  the  corresponding 
equations  are  also  given  in  the  figure. 


5.  SIMULATION  RESULTS 

A  two-dimensional  single  target  tracking 
problem  will  be  considered.  Two  target  dynamic 
models  will  be  assumed,  one  with  (nearly) 
constant  velocity  and  the  other  with  (nearly) 
constant  acceleration.  The  Markov  transition 
matrix  of  the  models  is  known  and  given.  The 
initial  target  state  estimate  is  given  and  the  initial 
probabilities  of  the  two  target  models  are 
assumed  equal. 

The  target  dynamic  models  with  discretization 
over  time  intervals  of  length  T  are 

*(*)  =  W*)]*(*-l) 

+  G[M(k)]v(k  -  1)  (38) 
where  for  model  1,  the  nearly  constant  velocity 
model,  the  state  is 


and 


F  = 


T 

1 

0 

0 


y  ]' 

o 

0 

T 

1 


(39) 


(40) 


G  = 


0 

0 

T2I2 

T 


(41) 


T  T2n 
T 
0 
0 

The  process  noise  v(k)  =  (u, ,  yv]'  representing 
the  acceleration  during  one  period  is  a  zero 
mean  Gaussian  white  noise  vector  with 
covariance 

[?.,  o 

'  o 

For  models  2  (with  acceleration),  the  state  is 


x  =  \x 


y  y  y\ 


(42) 


and 


F  = 


1 

0 

0 

0 

0 

Lo 


a  = 


T2/ 2 
T 
1 
0 
0 
0 

T2/ 2 
T 
1 

0 

0 

0 


0 

0 

0 

T 

1 

0 


0 

0 

0 

T2/ 2 
T 
1 


(43) 


0 

0 

0 

T2I2 

T 


(44) 


where  the  process  noise  v(k )  representing  here 
the  acceleration  increment  over  one  period  is  a 
zero  mean  Gaussian  white  noise  vector  with 
covariance 

[42,  0 

L  0  qly.’ 

Assuming  only  position  measurements  to  be 
available,  then,  for  node  / 


z'(k)  =  H'x(k)  +  w,(/fc) 

where 

r  1  0  0  0  0  ()' 

~  LO  0  0  1  0  0. 


(45) 


(46) 


and  w‘(k)  is  a  zero  mean  Gaussian  white  noise 
vector  with  covariance 


'  G  o ' 

.0  r\.' 

To  overcome  the  fact  that  one  has  different 
state  dimensions  the  lower  dimension  vector  was 
augmented  with  suitable  zero  components 
(which  then  have  mean  and  variance  zero)  to 
make  it  compatible  with  the  higher  dimension 
state. 

With  sampling  interval  T  -  Is,  the  true  target 
is  simulated  with  constant  velocity  for  the  first 
seven  scans,  then  switches  to  constant  acceler¬ 
ation  for  the  next  seven  scans,  and  finally  returns 
to  constant  velocity  for  another  seven  scans.  The 
initial  target  state  is  assumed  to  be  [100  m, 
30  m  s'1, 0,  100  m,  15ms~',0]  and  the  acceler¬ 
ation  is  assumed  to  be  5  and  -5  m  s~:  for  the  x 
and  y  coordinates,  respectively. 

The  variances  of  the  process  noise  are  taken  as 
<?i.x  =  Qi.y  —  0. 1  (m  s-2)2  for  model  1,  the  nearly 
constant  velocity  model,  and  q2t  =  = 

1.0  (ms-2)2  for  model  2,  the  nearly  constant 
acceleration  model.  The  detection  probabilities 
for  both  sensors  are  equal  to  0.67  and  the  false 
alarm  rates  are  0.0001  m~2.  The  standard 
deviations  of  the  measurement  errors  are 
assumed  to  be  V(10)m  for  both  x  and  v 
coordinates  of  the  two  sensors.  The  Markov 
transition  matrix  for  the  model  parameters  is 
assumed  to  be 


0.9  0.1 
.0.1  0.9.' 

The  initial  state  estimate  is  generated  randomly 
with  mean  the  same  as  the  true  target  state 
and  covariance  matrix  equal  to 

diag [100,  1,0.1,  UK),  1,0.1], 

Three  different  configurations  will  be  tested. 
First,  each  sensor  will  be  simulated  indepen¬ 
dently  using  the  MMPDA  algorithm  described  in 
3c<_ik)ii  3.  Second,  a  centralized  processing  with 
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Trajectory  window 


Fig.  3.  Tracking  results  with  sensor  1  only  (one  sample  run). 


measurements  from  both  sensors  will  be 
simulated  using  the  same  MMPDA  algorithm. 
Finally,  the  distributed  case  will  be  simulated.  In 
this  case,  the  two  nodes  will  communicate  every 
scan.*  At  each  scan,  each  node  will  process  its 
own  sensor  measurements  first,  then  send  the 
local  processed  results  to  the  fusion  node.  After 
receiving  the  information  from  both  local  nodes, 
the  fusion  node  will  use  the  fusion  algorithm 
derived  in  the  previous  section  to  construct  the 
global  estimates  and  send  the  results  back  to 
^ach  local  node. 

Simulations  were  carried  out  with  50  Monte 
Carlo  runs.  The  results  of  one  sample  run  are 
shown  in  Figs  3-5.  Figures  3  and  4  show  the 
estimated  and  true  trajectories  of  the  target  with 
sensors  1  and  2,  respectively.  Figure  5  shows  the 
results  for  the  distributed  case  where  the  two 
sensors  interchanged  their  processed  results.  As 
one  can  see,  the  single  sensor  processed  results 
have  poor  performance,  and  the  target  is  lost  in 
both  cases.  Figure  6  shows  the  probability 
trajectories  of  model  2  for  the  three  cases  as 
calculated  by  the  corresponding  state/model 
estimators.  As  can  be  seen  from  the  figures,  in 

*  This  Is  totally  equivalent  to  the  centralized  configuration 
but  has  the  advantages  of  redundancy  and  reliability  for  a 
DSN  system  This  configuration  can  also  be  used  with  a 
lower  communication  rate  (Chang  e<  at  ,  198fi) 


both  single  sensor  cases  the  algorithm  fails  to 
detect  clearly  the  switches  of  the  target  between 
two  models.  The  distributed  algorithm  not  only 
responds  faster  in  detecting  the  first  jump  of  the 
target  from  the  constant  velocity  mode  to  the 
constant  acceleration  mode,  but  also  successfully 
detects  the  end  of  the  acceleration.  The 
centralized  algorithm,  which  is  not  shown  in  the 
figures,  performs  exactly  the  same  as  the 
distributed  one. 

The  average  performances  for  the  three 
configurations  for  50  runs  are  given  in  Table  1. 
The  centralized  and  distributed  algorithms 
successfully  track  the  target  in  43  out  of  50  runs 
(“successful  tracking”  is  defined  when  the 
estimated  target  position  is  within  30  m  of  the 
true  target  position  for  the  last  three  scans). 
However,  out  of  50  runs,  sensor  1  alone  and 
sensor  2  alone  only  track  the  target  successfully 
in  27  and  30  runs,  respectively.  The  r.m.s. 
position  errors  for  those  successful  runs  are  also 
calculated.  Similarly,  the  centralized  and  distrib¬ 
uted  algorithms  perform  better  than  the  single 
sensor  configurations.  Note  that  the  quality  of 
the  estimation  using  two  sensors  in  terms  of 
mean  square  error  is  significantly  better  than 
using  a  single  sensor. 

The  centralized  case  yields  the  upper  bound  of 
the  performance  for  the  distributed  configur- 
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Fir;.  5.  Tracking  results  of  distributed  case  (one  sample  run) 
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Fig.  6.  Model  2  probability  trajectories. 


ation  when  the  nodes  communicate  every  scan. 
The  simulation  shows  that  the  results  of 
the  distributed  algorithm  are  the  same  as  in  the 
centralized  algorithm,  which  confirms  the  theor¬ 
etical  equivalence. 

6.  CONCLUSION 

A  recursive  estimation  algorithm  that  accounts 
for  the  uncertainties  of  both  measurement 
origins  and  system  models  in  a  distributed 
framework  has  been  derived.  The  distributed 
estimation  technique  has  been  adopted  together 
with  the  probabilistic  data  association  (PDA) 
filter  in  conjunction  with  the  interactive  multiple 
model  (IMM)  scheme.  The  resulting  algorithm 
can  be  applied  to  track  a  maneuvering  target  in  a 
cluttered  tnviionment  with  distributed  sensors. 
Simulation  results  show  the  expected  perform¬ 


ance  of  the  algorithm.  With  full  communication 
rate,  the  distributed  case  performs  exactly  the 
same  as  the  centralized  case,  which  confirms  the 
theoretical  equivalence,  but  has  the  advantages 
of  increased  reliability. 
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Substituting  these  into  (3.10),  and  using  (E-2),  we  obtain 

Urn  e,(()  =  0.  (3.13) 

The  estimation  property  (E-3),  the  uniform  boundedness  ol y(t)  and  u(t), 
ind  (2.5)  the  definition  of  f  *,  imply  that 

lim  e(t)  =  0. 

Substituting  this  into  (3.11)  and,  again,  using  (E-2)  we  obtain 

limei(t)  =  0.  (3.14) 

Since  is  a  stable  polynomial,  we  can  establish  ii)  by  substituting 

(3.13)  and  (3.14)  into  (3.12).  VVV 

Remark  3.1:  The  multirate  sampling  estimation  algorithm  in  general 

does  not  have  the  property  that  e(f)/[I  +  ||^(/  -  l)||2Jl/l  €  /j,  which  is 
required  in  the  stability  proof  of  conventional  adaptive  control  algorithms. 
However,  we  still  prove  the  stability  using  property  (E-3)  and  the  relation 
|e(/*)|  2:  (<?(/)(  for  £/</;• 

IV.  CONCLUSIONS 

In  this  note,  we  have  developed  a  multirate  sampling  adaptive  control 
algorithm  which  allows  a  fast  sampling  rate  of  feedback  control  to  be  used 
even  if  the  computation  of  parameter  estimate  and  controller  coefficient 
may  take  a  relatively  long  period  of  time. 

The  key  idea  to  achieve  this  is  to  record  the  plant  input  and  output  prior 
to  the  currently  obtained  estimate  and  use  them  to  compute  the  coming 
estimate  and  controller  coefficients.  Thus,  the  computation  is  not 
dependent  upon  the  inputs  and  outputs  appearing  during  the  updating 
process.  The  closed-loop  system  is  shown  to  be  stable. 

Remark  4. 1; 

i)  One  may  further  extend  the  algorithm  to  consider  tj  -  _  i  >  n  +  m 

+  d  =  ii.  In  this  case,  a  relation 

\eUj-t  +  /I  +  ir)(sCi  max  (e(<)l  +  Cj 

-  * 

(k  <  o°,  C,  <  oo,  Cj  <  co),  can  be  used,  and  the  algorithm  only  needs  to 
compute  e(r)  for  ty-i  £  t  <  +  rf  but  not  for  every  /  in  |  s:  t  <  tj. 

ii)  Instead  of  the  ARMA  model,  one  can  use  6-model  [8]  in  the 
algorithm,  which  retains  the  key  features  of  the  continuous-time  model 
and  allows  a  wide  bandwidth  MRAC  system  to  be  achieved. 

iii)  The  multirate  sampling  adaptive  control  is  presented  for  an  indirect 
MRAC  system.  However,  the  method  covers  a  wide  class  of  direct  and 
indirect  adaptive  control  algorithms  of  certainty  equivalence  type  such  as 
pole-assignment,  LQ -optimal,  etc. 

iv)  Various  methods  developed  for  improving  adaptive  control  system 
performance  are  applicable  to  the  presented  multirate  sampling  adaptive 
algorithm.  These  methods  include:  a)  various  modifies  -  ons  of  parameter 
estimator  for  improving  convergence  rate;  b)  noise  and  disturbance 
filtering  techniques;  c)  robustness  techniques  with  respect  to  disturbances 
and  unmodeled  dynamics,  such  as  deadzone,  normalization,  etc.;  d) 
internal  model  principle  for  deterministic  disturbance  rejection,  etc. 
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An  Adapiive  Dual  Controller  for  a  MIMO-ARMA 
System 

P.  MOOKERJEE  and  Y.  BAR-SHALOM 

Abstract— An  adaptive  dual  controller  Is  presented  here  for  a  multiin- 
pul  multiontput  ARMA  system.  The  plan!  bas  constant  but  unknown 
parameters.  The  cautious  controller  with  a  one-step  horizon  and  a  new 
dual  controller  with  a  two-step  horizon  are  examined.  In  many  Instances, 
the  myopic  cautious  controller  Is  seen  to  (urn  off  and  converges  very 
slowly.  The  dual  coatroller  modifies  (he  cautious  con(rol  design  by 
numerator  and  denominator  correction  terms  which  depend  upon  (he 
sensitivity  functions  of  (he  expected  future  cost  and  avoids  (he  turn-off 
and  slow  convergence.  Monte-Carlo  comparisons  based  on  parametric 
and  nonparamctric  statistical  analysis  indicate  (he  superiority  of  the  dual 
controller  over  the  cautious  controller. 

I.  Introduction 

Multiinput  multi output  systems  with  unknown  parameters  are  encoun¬ 
tered  in  many  practical  situations,  and  their  control  poses  a  great 
challenge  to  the  stochastic  control  theory.  It  is  not  possible  to  obtain  an 
optimal  solution  for  such  systems  because  of  the  dimensionality  involved 
in  the  stochastic  dynamic  programming  [6],  In  such  situations,  emphasis 
is  on  obtaining  a  suboptimal  solution  that  incorporates  the  intrinsic 
properties  of  the  optimal  solution.  For  stochastic  systems,  the  control  has 
in  general  a  dual  effect  p],  [1 IJ:  it  affects  the  system's  slate  as  well  as  the 
future  state  and/or  parameter  uncertainty.  Thus,  the  dual  controller  offers 
significant  improvement  potential  for  the  control  of  uncertain  linear 
plants.  In  multistage  problems  it  “probes'*  the  system  to  enhance  real¬ 
time  identification  of  the  system's  parameters  in  order  to  increase  the 
accuracy  of  the  subsequent  control  decisions  and  regulates  the  system  at 
the  same  time  {4],  [9], 

Two  classes  of  dual  controllers  exist  presently  [14J.  In  the  first  class 
[10],  [12],  [18],  the  control  minimizes  a  one-step  ahead  criterion 
augmented  by  a  second  terra  which  penalizes  for  poor  identification.  This 
approach  is  simple  but  often  requires  tuning  of  some  parameters.  The 
second  class  (developed  for  SISO  systems  in  [3],  [16],  [17])  used  the 
stochastic  dynamic  programming  equation  and  expands  the  future  cost 
about  a  nominal  trajectory.  Using  first-  and  second-order  Taylor  series 
expansions  of  the  expected  future  cost  about  a  nominal  trajectory,  dual 
controllers  for  MIMO  sucic  systems  are  developed  in  [5]  and  [14].  A 
second-order  Taylor  series  expansion  of  the  future  expected  cost  is 
performed  about  a  nominal  trajectory  and  a  dual  controller  based  on  • 
two-step  horizon  is  developed  in  this  note  for  a  MIMO  dynamic  (ARMA) 
model.  The  cautious  [14],  [16],  [18]  and  the  new  dual  controller  are 
applied  to  a  MIMO-ARMA  system.  Monte  Carlo  simulations,  based  on 
parametric  and  nonparametric  statistical  analysis,  indicate  that  the  dual 
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controller  prevents  the  turn-off  phenomenon  and  slow  convergence 
prevalent  with  a  cautious  solution. 

Section  □  gives  the  problem  formulation.  The  approximate  dual 
controller  with  a  two-step  horizon  for  the  MIMO  system  is  derived  in 
Section  m.  The  control  solution  is  obtained  by  approximating  the  solution 
of  the  stochastic  dynamic  programming  equation.  A  second-order  Taylor 
series  expansion  of  the  expected  future  cost  is  performed  about  a  nominal 
trajectory  and  this  leads  to  a  dual  control  solution  in  a  closed  form. 
Following  the  derivations  of  the  controller,  a  summary  of  the  algorithm  is 
given.  Section  IV  describes  the  simulation  of  the  plant  and  compares  the 
performances  of  the  cautious  and  the  dual  solutions.  Section  V  concludes 
the  note. 

II.  Problem  formulation 
The  MIMO  system  to  be  controlled  is  described  by 

y(k)= -Ay(k-l)  +  Bu(k-l)  +  e(k)  (1) 

where 

£[e<*)]  =0;  E[e(k)  =  Witj.  (2) 

Here  y(k)  is  the  output  of  the  plant,  u(k)  is  the  input  to  the  plant,  and 
e[k)  is  the  measurement  noise. 

The  parameter  matrices  A  and  B  are  unknown.  This  model  describes 
some  industrial  processes  like  an  ore  crushing  plant,  or  a  heat  exchanger 
(1].  The  unknown  elements  of  A  and  B  comprise  the  parameter  vector 
6(k)  whose  estimate  at  time  k  is  6(k)  with  covariance  matrix  P(k).  The 
parameter  vector  is  designated  as 

*(*)  6  [o.'lft.'lfljlhjl  •••  Klb.T  (3) 

where  n  is  the  dimension  of  the  output  vector  y(k)  and  a,' ,  b,'  are  the  / th 
row  of  the  matrices  A  and  B,  respectively.  Assuming  the  parameters  are 
time-invariant,  we  have 

6(k+l)~6(k).  (4) 

.  A  measurement  matrix  H(k)  is  defined  as 

H(k)  &  diag  [-y-'Wlu'W.  -y'(k)\u-(k),  ■  ■  •]  (5) 

where  H(k)  has  n  rows,  and  y'(k),  u'(k)  are  the  measurement  and 
control  vectors  transposed. 

With  these  definitions,  the  measurement  model  is 

y(k)  =  H(k-l)e{k-l)  +  e(k).  (6) 

The  performance  criterion  to  be  minimized  is  /(0),  i.e.,  the  conditional 
expected  value  of  the  cost  C(0)  from  step  0  to  N,  denoted  by 

y(0) =£•{0(0)1/*} 

where  Q(k)  is  the  diagonal  weighting  matrix,  /*  is  the  cumulated 
information  at  time  k,  and  y,  is  the  desired  output. 

in.  Dual  Control  with  a  Two-Step  Horizon 

First  the  controller  is  derived  and  then  a  summary  of  the  algorithm  is 
provided. 

A  dual  control  solution  with  a  two-step  horizon  is  obtained  by 
minimizing  (2.7)  with  respect  to  the  control  u( 0)  for  the  multidimensional 
plant  (2.t)-(2.4).  This  is  obtained  by  solving  the  general  equation  of 
stochastic  dynamic  programming  [3],  (7],  (8] 

7*(k)  =  min  £{C(k)  +  /*(k+  1  )| /* }  *  =  1,  ■••,1.0  (1) 

»<*) 

where  /*(k)  is  the  optimal  expected  cost  to  go  from  k  to  N,  C(k)  is  (he 


cost  to  go  from  k  to  N,  and  /*  is  the  cumulated  information  at  time  k  when 
the  control  u(k)  is  to  be  applied.  The  information  /*  is  the  set  of  all  past 
controls  until  time  k  -  I  and  outputs  until  time  k. 

Thus,  for  a  two-step  horizon  we  have 

-min  £[(>-(*+ l)-y,)'C(k){y(k+ l)-y,}+/;  |/»] 

(2) 

where  J*t ,  ttl  is  the  optimal  expected  cost  at  the  last  step  with  one -step 
horizon  and  is  obtained  by  minimization  of  Jt,  M,2,  ar|d  ./*,  2  is  the 

cost  to  go  from  k  +  1  to  k  +  2. 

The  cautious  control  at  Ar  +  1  with  one-step  horizon  is  given  by 

u(*+  l)  =  (£{£'Q(k+  1)£|/***}]  ” 1 

•  E[B'Q(.k  +  \  ){Ay(k  +  1)  +y, }  |f*  * ']  (3) 

The  cost  from  step  k  +  1  to  A  -t-  2  is 
7tt|.l.l  =  'fW+l)H' 

+  £[{/ty(k+  l)+y,)'Q(k+  l){Ay(k+l)  +  y,\ 

+  «'(*  +  l)B'Q(k  +  l)Bu(k  +  1)  -  2  {Ay(k+\)+y,)' 

■  Q(k+\)Bu(k+\)\lk"]  ■  (4) 

and  inserting  (3)  into  (4)  the  optimal  cost  at  the  last  step  is 

'*\i.*-2=«re(*+u»' 

+  £[{Ay(k+  l)+^}'Q(*-t-  IH^CAr-i- 
-£ri{Ay(fc+  l)  +  yr}‘Q(k+  1)B|/*'") 

•  (£{£'(>(*+ 1)B|/**'}]'' 

•  E[B'Q(k+l){Ay{k+l)+yr}\I^']'  (5) 

where  £■  { - 1 /*  ^ 1 }  is  the  conditional  expectation  given  the  available 
information  I** '. 

The  unknown  parameters  will  be  chosen  from  the  Gaussian  family  and 
thus  their  estimate  9(k  +  1)  and  associated  error  covariance  P(k  +  l)are 
the  sufficient  statistic.  The  parameter  vector  estimate  S(k  +  1)  and  the 
associated  covariance  matrix  P(k  +  1)  are  obtained  from  a  Kalman  filter 
according  to 

K(k+  1  )  =  P(k)H'(k)[H(k)P(k)H'(k)+  W]-'  (6) 

S(k+  1)  =  S(*)  +  A(k+  1  )[y(k+  l)-//(*)9(*)] 

=  Hk)  +  K(k+l)>’(k+l)  (7) 

P(k+  \)  =  P(k)-P(k)H\k)[H(k)P(k)H'(k)+  W)-'  H(k)P(k).  (8) 

Here  v(k  +  1)  is  the  innovation  of  the  process. 

From  (5)  it  is  clear  that  J*+ ,  2  is  a  nonlinear  function  of  the  estimated 

parameter  vector  S(k  +  I)  and  covariance  P(k  +  1).  But  the  estimated 
vector  9(k  +  1)  and  the  covariance  P(k  +  1)  are  not  known  until  the 
control  u(k)  is  applied. 

A  control  u(k)  with  a  two-step  horizon  can  be  obtained  from  (2)  if  a 
second-order  Taylor  series  expansion  of  t  +  2 's  performed  about  a 
suitable  nominal  trajectory.  Here  the  nominal  trajectory  is  defined  by 

1)  a  nominal  parameter  estimate  $(k+  1)  =  9(k) 

2)  a  nominal  control  tl(Ar) 

3)  a  nominal  covariance  P(k+  1)  obtained  by  using  d(k) 

4)  a  nominal  measurement  y(k  +  1)  obtained  by  using  il (k)  and 

$(k),  i.e..  y(k+l)  =  /}(k)S(k). 
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Expansion  of  (5)  about  this  nominal  trajectory  results  in 
i,  =}>  +  J',  (*+!)(>«:+  1) -y<*+  Di 

+  j  [>(*+  «)-/(*+  I)]'/,,(*+l)|j>(*  +  l)-?(k  +  1)1 
+  y;(Ar+l)[5(Jr+l)-J(/t)]+i[3«r+l)-«(Ar)l- 


with  the  superscript  here  denoting  the  matrix  element,  e,  the  t'th  Cartesian 
basis  vector,  and 


/>;;(*+  0 


&  3PHk+l) 
~  3u(k ) 


Pu„(k+l) 


3uJ(*) 


r,y« I,  •••./■ 


(70) 


•  y„(*+l)(3(*+l)-«(*)] 

+  tr  (y,,(*+  !){/>(*  +  l)-/5(Jt+  1)}1 
►here  7,  is  the  zeroth -order  term  and  the  cost  sensitivities  are 


(9) 


(10) 


(U) 

(12) 


*(<r+  0  =  |_a3,(*+i)aSy(*+i) 


(13) 


y,(<r+l)  =  3J —  .  (14) 

|_aP't(JI:+l)  J 

The  above  sensitivities  are  evaluated  at  d(k),  P(k  +  1),  and  y(k  +  1); 
utd  P<>(k  +  1)  is  the  yth  element  of  the  covariance  matrix  associated 
with  the  parameter  estimates  S,(k  +  1)  and  Sj(k  +  1). 

Under  the  Gaussian  ass„.nption  for  the  zero  mean  noise 


y(k+  l)-y(k  +  1)- 31 W.  V]  (15) 

►here  the  conditional  mean  is 

:hL^E{H(k)6(k)  +  e(k  + 1  -f?(*)9(*)|/*} 

=  [H(k)-/Hk)]d(k)  (16) 

ind  the  conditional  covariance  is 


'K=E[{y(<r+I)-jJ(<r+l)-,r}{y(Ir+l)->'(<r+l)-r<}'|/*l 

=  H(k)P(k)H'(k)+W.  (17) 

With  the  choice  of  the  nominal  path  as  defined  earlier  and  using  (6), 
(16),  and  (17),  the  conditional  expected  value  of  (9)  is 

+  iMV„(lfc+I)j.  +  i«r  V„(k+  !>n 
+  ^  tr  [J*(k  +  l){P(k)-P(k  +  1)}] 

+  tr  [Jr(.k+  I){P(E+  !)-/*(<:+  I)}J.  (18) 

The  above  expected  future  cost  (18)  is  a  function  of  the  nominal 
parameters  multiplied  by  appropriate  sensitivity  functions  J,{k  +  I), 
l„(k  +  1),  Jn(k  +  1),  and  J/fk  +  1).  These  sensitivities  introduce  the 
4ual  effect  into  (2)  which  is  then  used  to  yield  u(k).  It  must  also  be  noted 
hat  the  covariance  P(k  +  1)  is  nonlinear  in  u(k)  and  is  not  yet  known. 
Hence,  a  second-order  expansion  of  P(k  +  1)  is  proposed  about  a 
•ominal  control  U(k)  and  a  nominal  covariance  P(k  +  1)  in  order  to 
obtain  a  (suboptimai)  dual  solution  uD(k)  in  a  closed  form  from  (2). 
This  expansion  is  performed  as  follows: 

/>(*+  |)  =  P(*+  1)+  Yj  e.</  ^PIJ.(k+  !)(</(*)- <7(*)1 

+  ^  [«(*)- a(k)]-pl(k+  l)[u(4r)-d(4r))]  (19) 


evaluated  at  P(k  +  1)  and  ii(Ar)  and  r  the  number  of  unknown  parameters. 

Now  a  (suboptimai)  dual  solution  uD(k)  with  a  two-step  horizon  can  be 
obtained  from  (2)  using  (l8)-(20)  and  is  given  in  closed  form  by 

uD(k)  =  {E{B-Q(k)B\l‘\  +  F)-'[E{B'Q{k)(Aylk)+y,)\P}+f]  (21) 

where  the  elements  of  the  matr:  i  F  and  those  of  the  vector / are  given  by 


rU- 


u  j^y. 


(k+!)--y«(A:+l) 


+  2U 


+  2tf 


)  32P(k  + 1)  ~| 
j  3u,(k)3uj(k)  J 


J„ik+  1) 


V  3u,(k) 


i.  j  =  1 ,  •  ■  • ,  m 


(22) 


and 


-l.  ig^] 


1„ 


(23) 


and  m  is  the  dimension  of  the  control  vector,  ir,  is  the  nh  element  of  the 
control  vector. 

It  is  clear  from  (21)  that  this  approximate  dual  solution  uD(k)  is  a 
modification  of  the  cautious  solution  by  the  cost  sensitivity  terms.  The 
cautious  solution  is  (21)  with  F  -  Q  and  /  =  0.  These  account  for  the  dual 
effect.  The  implementation  of  this  second-order  dual  solution  is  per¬ 
formed  by  the  method  described  below. 

Algorithm  Summary: 

1)  Compute  the  sensitivity  functions  Jn(k  +  1),  J?(k  +  I),  J,(k  + 
1),  J„(k  +  1)  for  (18)  with  S(k  +  1)  =  6(k)  and  the  nominal  values 
tl(k),  P(k  +  I),  f(k  +  1)  defining  the  nominal  path. 

2)  Search  on  (2)  with  (18)  [with  the  sensitivity  functions  computed 
above,  starting  with  first  nominal  values  Q(k),  P(k  +  1)]  over  u(k)  to 
obtain  an  improved  nominal  for  which  f*  ktl  is  lower.  This  search  is 
done  by  selecting  a  first  coarse  grid.  A  grid  search  is  necessary  to  avoid 
locking  in  on  a  local  minimum.  Then  another  grid  is  chosen  about  the 
latter  control  over  a  narrower  interval  and  from  a  second  search  u‘(k)  is 
obtained. 

3)  Using  u'(k)  compute  the  covariance  sensitivities  P,(k  +  1 ) ,  (k 
+  1);  together  with  the  previously  computed  cost  sensitivities  yw(k  +  1), 
JP(k  +  1),  J„(k  +  1),  J,(k  +  1)  obtain  F,  f  defined  in  (22),  (23). 
Finally,  the  control  to  be  applied,  uD{k),  is  calculated  from  its  explicit 
expression  (21). 

The  iteration  described  in  step  2)  above  is  carried  out  to  obtain  better 
covariance  sensitivities.  The  control  Uo(k)  could  have  been  obtained 
directly  from  (21)  by  skipping  step  2)  above;  however,  as  indicated  in  [1 3] 
and  1 14],  this  results  in  unsatisfactory  performance.  With  this  iteration  of 
step  2),  the  “improved”  sensitivities  yield  good  performance  as  shown  in 
the  next  section. 
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rv.  Simulation  Results 

Performance  is  evaluaied  from  500  Monle  Carlo  runs  for  the  following 
controllers: 

1)  heuristic  certainty  equivalence  [3]  (with  a  one-step  horizon); 

2)  one-step  ahead  cautious  controller;  and 

3)  dual  controller  based  upon  sensitivity  functions  (with  a  two-step 
horizon)  derived  in  Section  ni. 

The  plant  equations  for  a  two-input  two-output  system  are 

y,(<r+l)=  -a„y,(k)-al2y2(k)  +  b„ut(k)  +  b,2u2(k)  +  e,(k  +  I)  (I) 
yi(k+  I)=  -a2lyt(k)-a22y2(k)  +  b2iu,(k)  +  bnu2(k)  +  e2(k  +  1)  (2) 
where 

C{c{k)e'(J)\  =  W/i»/  =  diag  (  W,,  W2)- 

H',  =  7.52J;  If,  =  43s.  (3) 

The  true  values  of  the  parameters  are 


Cn  =  0.8 

-74.84 

© 

II 

<5" 

b,2  = 

-51.04 

a2<  =  0.2 

53.31 

a„  =  0.75 

6„  = 

-82.56. 

(4) 

Only  the  gain  parameters  (B  matrix)  are  considered  unknown  for 
testing  the  dual  effect  and  their  initial  estimates  were  generated  as  31  (by, 
bV),  i,j  =1,2.  This  choice  of  system  was  motivated  by  the  helicopter 
vibration  study  [13], 

A  large  initial  uncertainty  is  chosen  in  the  parameter  estimates  in  order 
to  test  the  learning  capabilities  of  the  various  adaptive  algorithms.  The 


cost  weighting  matrices  are 

Q(fc)  =  diag(q„  q,);  q,  *  1.0,  q,=  1.0.  (5) 

The  desired  response  is 

3V  =  [  - 1 8  80] ' .  (6) 

For  the  model  chosen  (1)— <6)  the  optimal  control  solution  in  order  to 
reach  a  steady-state  value  of  y,  in  (6)  is 

u*  =  1.0,  uj=  -1.0.  (7) 

In  terms  of  the  notation  of  (1)  and  (2) 

9(k)  &  fa„  al2  £„(*>  £„(*)  <z„  aa  6u(k)  6n(k)Y  (8) 


TABLE  I 

average  Costs  for  the  Three  algorithms  in  the  simulation 

WITH  A  LIMITER  (|u,|  S  2.0,  | if, |  s  2.0)  (500  MONTE  CARLO  RUNS) 

the  superior  Rate  of  adaptation  of  the  dual  algorithm 
is  demonstrated  here 


Tlae 

Step 

HCE 

Cautipus 

Du*l 

k 

1 

k 

k 

& 

1-1 

53 

1-1 

h 

i-i 

I 

14851 

14851 

3623 

3623 

6944 

6944 

2 

6241 

21092 

3961 

7584 

6722 

13666 

3 

3578 

24670 

3246 

10830 

*230 

17896 

4 

1616 

26286 

2836 

13666 

1866 

19762 

5 

1354 

27640 

2505 

16171 

1492 

21254 

' 

807 

28447 

2154 

1632S 

953 

22207 

7 

593 

29040 

1921 

20246 

700 

22907 

8 

462 

29502 

1670 

21916 

582 

23489 

9 

397 

29899 

1623 

23539 

535 

24024 

10 

347 

30246 

1327 

24866 

385 

24409 

40 

77 

34444 

281 

43810 

69 

29178 

table  11 

statistical  Significance  test  for  comparisons  of  the  Cautious 

AND  THE  DUAL  ALGORITHM  IN  THE  SIMULATION  WITH  A  LIMITER 
(| if ,  |  eS  2.0.  |u,|  s  2.0)  (500  MONTE  CARLO  RUNS) 


ru« 

Taat 

Estimated 

Stap 

Statistic 

Improvement 

k 

1 

-e.i 

-91 

2 

-5. 3 

-69 

3 

-2.2 

-30 

4 

3.5 

34 

5 

3.3 

40 

6 

6.0 

56 

7 

6.3 

64 

8 

6.5 

65 

9 

6.5 

67 

10 

5.7 

71 

11 

6.3 

76 

12 

5.6 

70 

13 

5.9 

82 

14 

5.2 

62 

15 

5.5 

79 

16 

4.9 

70 

17 

4.5 

78 

16 

4.4 

74 

19 

4.4 

76 

20 

4.3 

76 

and 


H{k)  & 


-y,(.k) 

0 


-yiW 

o 


u,(k)  u2(k)  0  0  0  0 

0  0  -y,(k)  -y2(k)  u,(k)  u2(k) 


(9) 


The  controllers  are  implemented  with  a  sliding  horizon  for  a  total  of  40 
time  steps.  The  evaluation  criterion  is 

Ci-Wt+D-hl'efW+ll-h)-  (10) 

A.  Analysis  of  the  Monte  Carlo  Average  Costs 

Comparisons  are  made  between  the  performances  of  the  cautious  and 
the  dual  algorithm  on  the  system  and  a  statistical  significance  analysis  is 
done  using  the  normal  theory  approach  (i.e.,  it  is  assumed  that  the  central 
limit  theorem  holds  for  the  sample  mean  from  a  large  number  of  runs) 
[14).  Tables  I-IV  contain  the  results  of  the  simulation  runs.  Table  I 
compares  the  average  cost  C2  over  500  Monte  Carlo  runs  for  the  first  40 
time  steps  for  .HCE,  cautious  and  the  dual  algorithms,  with  a  control 
limiter  |ir,|  s  2,  /  =  1,  2. 

Clearly  it  is  seen  that  the  cumulative  average  cost  is  the  lowest  for  the 
dual  controller.  The  HCE  incurs  an  excessive  penalty  in  time  step  1) 
because  of  lack  of  caution.  The  cautious  controller  is  overly  cautious  and 
exhibits  slow  convergence.  However,  the  dual  controller  incurs  less 
penalty  in  time  step  1)  than  the  HCE  and  makes  a  judicious  choice  of 


caution  and  probing  to  leam  the  parameters  fast.  Fig.  1  compares  the 
performances  of  the  three  algorithms  for  500  Monte  Carlo  runs.  Both 
Table  I  and  Fig.  1  demonstrate  the  superior  rate  of  adaptation  of  the 
dual  algorithm. 

Table  II  provides  a  statistical  significance  test  and  shows  the  improved 
performances  of  the  dual  solution  from  time  step  4)  onwards  with  at  least 
98  percent  confidence. 

Table  HI  indicates  the  percentage  of  runs  where  the  cost  exceeds  2000 
for  the  two  algorithms.  This  threshold  of  2000  is  selected  from  a  sample 
distribution  study  of  the  cost  at  each  time  step.  Table  IV  shows  the 
percentile  test  [14],  [15]  comparing  the  cautious  and  the  dual  solution. 
They  clearly  indicate  from  time  step  4)  onwards  the  light  tailed  nature  of 
the  distribution  of  the  cost  yielded  by  the  new  dual  control  algorithm. 

B.  Individual  Time  History  Runs 

Analysis  of  the  Monte  Carlo  average  cost  indicates  the  improvement 
offered  by  the  dual  solution;  it  provides  no  information  sbout  the  cautious 
control's  turning-off  phenomenon  [16],  [18],  Hence,  •  careful  investiga¬ 
tion  of  the  individual  runs  is  required  to  examine  these  occurrences. 
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TABLE  111 

COMPARISON  OF  THE  TAILS  USING  THE  CAUTIOUS  AND  THE  DUAL 
ALGORITHMS  IN  THE  SIMULATION  WITH  A  UM1TER 
(|U||  £  2.0,  |uj|  £  2.0)  (500  MONTE  CARLO  RUNS) 


T1m 

Stap 

P«rc«ntag*  of  run* 
which  axcaad  2000 

k 

Cautloua 

Dual 

1 

86 

76 

2 

60 

52 

3 

43 

40 

4 

33 

25 

5 

•  31 

17 

6 

22 

10 

7 

22 

8 

6 

19 

7 

9 

16 

3 

to 

12 

2 

u 

12 

1.2 

12 

10 

1.4 

13 

11 

1.4 

14 

7 

1 

15 

8 

0.4 

16 

6 

0.4 

17 

6 

0.2 

16 

6 

0.4 

19 

5 

0.4 

20 

5 

0.2 

CAUTIOUS  AND  DUAL 


Fig.  2.  Time  history  of  output  1  using  the  cautious  and  the  dual  algorithms 
for  run  90  (500  Monte  Carlo  runs;  |u,|  £  2.0;  |u2|  £  2.0). 


f  TABLE  IV 

PFP.CENTLE  TEST  FOR  COMPARISONS  OF  THE  CAUTIOUS  AND  THE  DUAL 
I  ALGORITHMS  IN  THE  SIMULATION  WITH  A  LIMITER 

i  (|t/,|  s  2.0.  | Ui |  <;  2.0)  (500  Monte  Carlo  runs) 


Tits* 

Stap 

k 

X 1  caat  acaclatlca 
*c  K.0 

1 

2 

3 

4 

10 

5 

19 

6 

23 

7 

32 

8 

35 

9 

57 

10 

37 

11 

•40 

12 

40 

13 

40 

14 

16 

15 

32 

16 

11 

17 

16 

18 

16 

19 

18 

20 

25 

Time  St«p 

Fig.  3.  Time  history  of  output  2  using  the  cautious  and  the  dual  algorithms 
for  run  90  (500  Monte  Carlo  runs;  |ui(  £  2.0;  |«j|  £  2.0). 


CAUTIOUS  AND  DUAL 


CAUTIOUS.  DUAL  AND  HCE 


Fig.  1.  Time  history  of  the  average  cost  using  the  heuristic  certainty 
equivalence,  cautious,  and  the  dual  controllers.  (500  Monte  Carlo  runs; 
|«r,|  £  2.0,  |u2|  £  2.0.)  The  superior  rate  of  adaptation  of  the  dual 
algorithm  is  demonstrated  here. 


Fig.  4.  Time  history  of  control  I  using  the  cautious  and  the  dual  algorithms 
for  run  90  (500  Monte  Carlo  runs;  |u,|  s  2.0;  |tr2|  s  2.0). 
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Fig.  5.  Time  history  of  control  2  using  the  cautious  and  the  dual  algorithms 
for  run  90  (500  Monte  Carlo  runs;  (ut|  s  2.0;  |tr2|  s  2.0). 

The  turn-off  phenomenon  is  observed  in  many  runs  among  the  500 
Monte  Carlo  simulations  while  using  the  cautious  controller;  run  90  is  a 
typical  example  of  it.  Both  components  are  almost  off  between  time  steps 
0  and  20  during  which  the  dual  controller  already  identified  the 
parameters  and  reached  the  desired  trajectory.  Figs.  2-5  portray  this 
result. 

V.  conclusions 

A  new  adaptive  dual  control  solution  with  a  two-step  sliding  horizon 
has  been  developed  for  an  ARMA-MIMO  system.  The  control  law  is 
derived  by  solving  the  stochastic  dynamic  programming  equation.  This 
solution  utilizes  the  dual  effect  by  performing  a  second-order  Taylor 
sejies  expansion  of  the  expected  future  cost  and  does  not  need  any  tuning 
for  any  of  the  runs  in  the  example.  It  modifies  the  cautious  solution  by 
explicit  numerator  and  denominator  correction  terms.  The  controller  in  its 
present  form  is  the  first  of  its  kind  in  a  closed  form  for  a  system  with 
unknown  parameters.  The  controller  is  tested  on  a  MIMO  system  in  a 
systematic  Monte  Carlo  fashion.  Conclusions  are  based  on  500  Monte 
Carlo  nins..  Analysis  of  the  simulation  runs  has  shown  that  this  new  dual 
control  solution  applied  to  a  multiinpui  multioutput  model  improves 


over  the  cautious  controller.  The  key  improvement  is  in  the  avoiding  of 
situations  like  turn-off  and  stow  convergences,  typical  of  the  cautious 
solution. 
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ABSTRACT 

Tha  ravaralon  in  tiaa  of  a  atochaatic  dlffaranca 
aquation  In  a  hybrid  apaca,  with  a  Markovian 
solution,  ia  praaantad.  Tha  ravaralon  la  obtainad 
by  a  martingale  approach,  which  pravioualy  lad  to 
ravaraa  tiaa  forma  for  atochaatic  equations  with 
Cauaa-Harkov  or  diffuaion  aolutiona.  Tha  ravaraa 
.tiaa  aquationa  follow  froa  a  particular 
non-canonical  aartinqala  decomposition,  vhila  tha 
ravaraa  tiaa  aquationa  for  Cauaa-Karkov  and 
diffuaion  aolutiona  followad  froa  tha  canonical 
aartinqala  dacoapos ltion.  Tha  naad  for  thia 
non-canonical  daconpoaition  ataaa  froa  tha  hybrid 
•tata  apaca  aituation.  Horaovar,  tha  non-Cauasian 
diacrata  tiaa  aituation  laada  to  ravaraa  tiaa 
aquationa  that  incorporata  a  Bayaaian  aatiaation 
atap. 

1.  INTRODUCTION 

This  papar  adressaa  tha  problaa  of  tlaa-ravaraion 
of  a  hybrid  atata  Markov  procaaa  which  ia  qivan  aa 
tha  aolution  of  a  atochaatic  dlffaranca  aquation. 
Tha  daairad  raault  la  a  alailar  aquation  but 
runninq  in  ravaraa-tlna  diraction  vhila  havinq  a 
aolution  that  ia  raapactivaly  pathwiaa  and  in 
probability  law  aquivalant  to  tha  aolution  of  tha 
forward  aquation. 

Tha  motivation  to  atudy  this  problaa  ataaa  froa 
two  dlffarant  kinds  of  application.  Tha  first  is 
to  approach  tha  aolution  of  a  nonlinaar  saoothlnq 
problaa  by  a  aarqlng  of  tha  astiaatas  of  two 
nonlinaar  f lltars t  ona  fiitar  aatchaa  tha  oriqinal 
modal  and  is  appllad  in  tha  usual  tiaa  diraction 
vhila  tha  othar  fiitar  aatchaa  tha  tiaa-rsvarsad 
aodal  and  is  appllad  in  tha  ravarsa-tiaa 
diraction.  Tha  sacond  application  is  tha 
dataralnatlon  of  a  rata  distortion  thaory  lowar 
bound  for  a  diacrata-tlaa  nonlinaar  filtarinq 
problaa  by  tha  aathod  of  Caldos.  This  aathod  la 
baaad  on  Bucy'a  raprasantatlon  formula  and 
raquiras  a  Monta  Carlo  simulation  In  ravarsa-tiaa 
diraction  of  aodal  aatchlnq  trajectories,  startlnq 
from  a  praapaciflad  and  point  (Caldos,  19*1; 
Washburn  at  al.,  19*5).  Tor  both  of  thasa  two 
applications  it  is  nacassary  to  hava  a 
tiaa-ravarsad  dlffaranca  aquation  for  which  tha 
Markovian  solutions  ara  In  probability  law 
aquivalant  to  tha  oriqinal  solution. 
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Our  problaa  falls  in  tha  cataqory  of  how  to 
ravaraa  a  Markov  process  In  tiaa.  Tha  Markov 
proparty  laplias  that  tha  past  and  tha  future  ara 
lndapandant  undar  tha  condition  that  tha  praaant 
atata  is  known  (Wentzell,  19*1) .  This  lnvarlanca 
with  raapact  to  tha  tiaa  diraction  la  tha  kay 
proparty  uaad  in  tiaa-ravarslon  atudlaa.  Thara  ara 
two  typas  of  atudlaa  that  daal  with  thia  problaa; 
a  classical  typa  and  a  sy staas-typa .  Tha  classical 
typa  of  atudy  aaauaas  that  tha  transition  aaaaura 
or  tha  qanarator  of  a  Markov  procaaa  ia  qivan  and 
than  trlas  to  charactarita  tha  transition  aaaaura 
in  ravarsa-tiaa  diraction  (Naqasawa,  1964;  Kunlta 
and  Watanaba,  1966;  Chunq  and  Walsh,  1969;  Aztma, 
1971;  Kasaqava,  1976;  Oynkin,  197*;  Williaas, 

1979) . 

Tha  ayataaa-typa  of  atudy  aaauaas  that  a ’ 
stochastic  aquation  with  a  Markovian  solution  ia 
qivan  for  which  it  trlas  to  charactariza  tha 
tiaa-ravarsad  aquation.  Tha  first  tiaa-ravaraad 
aquations  wars  obtainad  by  orthoqonality 
arquaanta,  for  tha  linear  Causslan  situation 
(Ljunq  and  Kailath,  1976;  Lainiotis,  1976) .  For 
qanaral  diffusions,  it  has  alraady  baan  pointed 
out  by  Stratonovich  (1960)  how  to  obtain  tha 
ravaraad-tiaa  aquations  by  actually  followinq  tha 
classical  approach:  froa  a  stochastic  aquation  via 
tha  qanarator  and  tha  tiaa-ravarsad  qanarator  back 
to  tiaa-ravarsad  aquations.  A  truly  ayataas-typa 
of  study  has  baan  started  by  Verqhese  and  Kailath 
(1979),  by  showinq  how  for  a  linaar  Caua-ia.. 
system  a  more  direct  aartinqala  approach  laada  in 
a  slaplar  way  to  tiaa-ravarsad  aquations. 

Moreover,  by  this  approach  it  was  possible  to 
obtain  a  ravaraad-tiaa  aquation  with  a  pathwiaa 
aquivalant  solution.  Early  elaborations  of  these 
ldaas  lad,  alonq  different  routes,  to 
time-ravsrsed  aquations  with  pathwiaa  aquivalant 
aolutiona  (Anderson,  19*2;  Casta non,  19*2; 

Fardoux,  19*3).  Durinq  subsequent  studies,  quits 
larqe  classes  of  stochastic  differential  aquations 
and  thair  ravaraad-tiaa  aquations  hava  bean 
identified  (Elliott  and  Anderson,  19*5;  Fardoux, 
19(5;  Elliott,  1916a,  1916b;  Haussmann  and 
Fardoux,  19(6;  Pardoux,  19*6).  Recently  these 
results  hava  bean  extended  by  usinq  tha  Clrsanov 
transformation  of  Brovnlan  motion  (Picard,  19*6; 
Frottar,  19*7) .  Obviously,  this  Cirsanov  approach 
can  not  be  applied  to  discontinuous  or 
discrata^time  processes. 

To  qlva  an  idea  of  why  thara  is  an  additional 
problaa  in  usinq  a  aartinqala  approach  to  the 
raversiSh  of  an  aquation  with  a  discontinuous 
solution,  £e  qiva  a  brief  outline  of  tha  approach. 
Tha  aartinqala  approach  rouqhly  consists  of 
checkinq  if  tha  time-reversed  driving  noise 


sequence  can  ba  decomposed  In  a  suitable 
t«v*ca«-tlu  martingale  part  and  lta  complement 
and  next,  If  auch  a  decomposition  axlata  (Jecod 
and  Shlryaev,  ltt?;  Jecod  and  Prottsr,  ItaafT 
selecting  auch  a  decomposition.  Tha  final  atap  la 
to  characteclz#  both  tha  martingala  part  and  lta 
coaplaaant.  Xn  contraat  with  a  contlnuoua  procaaa 
auch  a  dacoapoaltlon  la  not  unlqua  for  a 
dlacontlnuoua  procaaa  (aaa  for  axampla,  Jacod  and 
Shlryaev,  19*7).  Thla  makes  tha  aalactlon  of“» 
aultabla  martingale  dacoapoaltlon  far  froa  trivial 
In  tha  hybrid  atata  apaca  situation,  bacauaa  a 
laaa  g&id  choice  yieade  unnacaaaarily  coaplicatad 
ravaraa-tlaa  equations.  Thla  complication  is 
prasantly  unaolvad,  naithar  In  continuous-tiaa  nor 
in  dlscrata-tiaa.  It  will  ba  aolvad  in  tha  saqual 
for  quits  ganaral  diffaranca  aquations  In  a  hybrid 
apaca.  With  that  rasult  wa  subaaquantly  ravaraa 
tha  conaidarad  diffaranca  aquation  In  time. 

Tha  papar  la  organized  as  follows.  In  aaction  2  wa 
dafina  tha  hybrid  atata  stochastic  diffaranca 
aquaClon  that  will  ba  conaidarad  and  shortly 
compara  its  tima-ravarslon  with  tha  time-reversion 
of  a  linaar  Cauaslan  aquation.  In  aaction  1  wa 
spaclfy  tha  tima-ravarslon  raqulrananta.  Next,  in 
sactions  4  and  S  wa  conaldar,  raapactivaly ,  tha 
pathwiaa  tima-ravarslon  and  tha  in  probability  law 
aqulvalant  tine-reversion.  In  aaction  4  wa  discuss 
tha  results  obtained. 

2.  THE  STOCHASTIC  DIFFERENCE  EQUATION  CONSIDERED 

The  stochastic  difference  equation  we  consider  In 
tha  sequel  is  tha  following  system,  on  an 
appropriate  stochastic  basis  and  a  discrete  time 
interval  [O.T]  *  Hx(0,TJ,  T<-, 

*t+l  "  «<*t+l«*t'xt.wt>  >  (!•»> 

«t+l  "  b(St,vt),  (l.b) 

yt  -  c(»t,xt,wt,ut) ,  (l.c) 

where  (vt),  (ut)  and  (vt)  are  i.i.d.  standard 
Gaussian  sequences  of  dimension  p,  q  and  1 
respectively,  the  initial  distribution  of  (x0,«0) 
has  tha  density  mass  function  p  ,  and 

*0<  *0 

(wt,vt,ut)  is  independent  of  (x0,*o>'  Further  xt, 
*t  and  yt  have  respectively  Rn-,  X-  and  R*-valued 
realizations  (with  K  a  countable  sat),  while  a,  b 
and  c  arm  measurable  mappings  of  appropriate 
dimensions  such  that  system  (1)  has  a  unique 
solution  for  each  initial  (Xg,*0)  with 

P„  .  (Xo.SoJ'O-  The  mappings  a,  b  and  c  are 
*0'®0 

time-invariant  for  notational  simplicity  only. 

The  second  order  dependence  of  (l.a)  on  («t)  is 
quite  uncommon  (Blom,  19*5).  Obviously,  (l.a) 
reduces  to  the  more  common  situation  of  first 
order  dependence,  only  If  a(«t«.i,*t>xt>vt)  is 
invariant  w.r.t.  either  «t  or  #t+1.  The 
interpretation  of  (l.a)  as  an  equation  with  a 
second  order  dependence  on  («t)  suggests  tha 
substitution  of  tt+i« (*t+l«*t>  ln  (*•»)•  °h  doing 
this  (l.a)  reduces  to  tha  more  common  equation, 
and  it  follows  immediately  that  (fct)  and  (f.t,xt) 
are  Martov  processes.  However,  as  the  state  specs 


of  4t  is  significantly  larger  than  the  atata  space 
of  at,  this  is  a  rather  brute  force  transformation 
of  (l.a).  A  mors  elegant  transformation  of  (l.a) 
to  the  more  common  equation  consists  of 
substituting  (l.b)  In  (l.a),  which  yields  an 
aquation  of  the  following  form, 

*t+l  *  •,(«t<*t'wt<vt)- 

Instead  of  a  state  space  expansion,  thara  appears 
an  additional  noise  term,  vt.  From  the  latter 
representation,  it  follows  immediately  that  the 
processes  («t,xc)  and  !*-)  •<-a  Markov  processes. 
The  latter  transformation  clearly  shows  that  (l.a) 
is  indeed  sore  general  than  tha  more  commonly 
studied  equation  with  first  order  dependence  of 
(at).  With  the  study  of  this  more  general 
aquation,  we  also  anticipate  tha  time-reversion 
results  obtained.  In  the  sequel  it  will  turn  out 
that  a  reverse-time  equation  of  (l.a)  has,  in 
general,  a  second  order  dependence  on  the 
time-reversed  (»t),  even  when  a(»t4.^,«t,xt,wt)  is 
•^-invariant.  In  view  of  this,  it  is  natural  to 
study  the  above  more  general  form. 

In  the  sequel  we  consider  the  time-reversion  of 
system  (1)  under  the  following  assumptions: 

lul 

a(#,s,.,w)  has  an  inverse  a*:KJxRnxRP«Rn,  such 
that  for  any  (« , s , w) £H2xrP, 

a*  (« ,s,a(»  ,k,x,wj',w)-xi  all  x£Rn.  (2) 

b(.,v)  has  an  inverse  b*:XxR-K,  such  that  for  any 
v€R, 

b*(b(s,v) ,v)-«;  all  »EK.  (3) 

Assumptions  A. 1  and  A. 2  suggest  to  transform 
(l.a.b.c)  to  the  following  time-reversed  model, 

t *t“*  (•t*l**t,xt+l'wtJ < 

(*t»l«vt) < 
yt“c(*t'*t>vt'ut)  • 

Because  (vt,vt)  and  the  future  (-  reverse-time 
past),  »t+1  -  »((y,, *«,»,);  S£(t+1,T)),  are 
dependent,  this  is  not  the  time-reversed  system  we 
should  look  for.  Unfortunately,  it  is  not  clear 
how  to  continue  from  here.  To  develop  some 
insight,  we  take  a  quick  look  at  the 
time-reversion  of  a  linear  Causslan  system. 


Consider  the  following  linear  Gaussian  system 
xt+l  *  **t  *  Bvf 

Assumption  A.  1  implies  that  A  is  invertible,  h-/ 
which 

xt  “  Cxt+1  "  ^tl- 

Obviously  W(  end  the  future  ere  dependent, 

which  requires  a  mart Inga Is  decomposition  of  wt. 
In  this  linear  Causslan  case  the  canonical 
martingale  decomposition  is  the  appropriate  one. 
It  consists  of  decomposing  vt  in  its  reverse-time 
predictable  part,  E(vt|»t+l),  and  its  complement 
w*» ; 

,  wt  -  ElWtISt*!)  ♦  v*t. 

The  problem  is  now  to  write  the  predictable  part 
as  a  function  of  xt*j  (if  possible)  end  to 


characterize  the  covariance  o f  v*t.  As  pointed  out 
by  Vergheee  and  Ka  1 lath  (1979)  it  follows  readily 
from  orthogonality  arguments  that 

e(wtl,telf  “  Elwtlxt«-ll< 

while  the  fundamental  formula  for  LLSE  estimation 
yields 

E«“tl*t+1»  “  *TR~1(tel)>u*1. 

Cov(w*t|  -  I-  8TR*l(t»l)B,  >. 

where  R(t+1)  is  the  covariance  of 
By  a  straightforward  substitution  of  these  remits 
we  obtain 

xt  -  A-»  (xt+1  -  B  BTR-^t+lJXt*!  -  Bw*t), 
which  yields  the  desired  reverse-tine  system: 

*t  -  A-l  (*t+1  -  B  BTR-1(t+l)*t4.1  -  BSt). 

The  orthogonality  arguments  and  the  LLSE 
estimation  step,  used  in  tb'.  above  procedure, 
prevent  a  straightforward  extension  of  that 
procedure  to  equation  (1).  In  the  sequel  we 
replace  the  orthogonality  arguments  and  the  LLSE 
estimation  step  respectively  by  Karkov  duality 
arguments  and  a  Bayesian  estimation  stap.  Besides 
this,  we  have  to  select  an  appropriate  martingale 
decomposition.  Following  the  linear  Cauasian  case, 
the  canonical  martingale  decomposition  seems  a 
good  candidate: 

(wt>vt)“<wt*»vt*)*E<  (wt<vtH»t+ll ' 
Unfortunately,  this  decomposition  leads  to  very 
complicated  elaborations  of  the  Bayesian 
estimation  stap.  To  avoid  these  complications,  we 
use  in  this  paper  the  following  decomposition: 

(«t**vt*)  "  <wt<vt)  “  («t*vtJ  * 
with:  vt  a  E(vt|»t+1)  and 
wt  s  E(vt|»t+1,vt). 

The  main  step,  that  must  be  carried  out,  is  to 
prove  that  the  latter  is  s  martingale 
decomposition,  and  to  elaborate  on  the  Bayesian 
estimation  step.  For  the  presentation  of  these 
results  a  constructive  approach  is  taken,  starting 
with  a  precise  description  of  the  time-reversio- 
objectives. 


*t 

»t 

?t 

where 


*(t«>tel«*t»*tel >°0  • 

(4. a) 

S(t.*t+l.*t41.*t> » 

<4.b) 

C (t, Bt+l# 0t* Et+1 ' Et * ®t* ut^  * 

(4.C) 

t,  B  and  6  are  deterministic  mappings  of 


appropriate  dimensions  and  *  n°lB* 

sequence  to  be  specified.  For  a  better 
understanding  of  (4)  notice  that  the  substitutions 
of  (4. a)  in  (4.c)  and  of  (4.b)  in  (4.a,c) 
transform  (4)  to  a  reversa-time  system  of  the  more 
common  form: 


it  “  ^^i^t+l'^Wl'^t'^t)' 

*t  “  B(t**t+l'*t+l*9t> > 

?t  “  e<t,*t>l,*t+l«Ot',t'ut)  >  t€(0,T-l]. 

To  be  a  useful  reverse-time  system,  (St,9t)  should, 
as  much  as  possible,  be  independent  of  the  future 
(*■  reve-sed-time  past)  information  field 

ffc+1  ■  • ( (?s* ES» ^s'^s'^s'^s) 7  s£(t+l,T]). 

A  minimal  requirement  is  then,  that  the 

conditional  expectation  of  (0t,9t),  given 

should  be  zero.  Because  9t  is  a  decreasing 
saquance  of  sigma  algabras,  the  latter  can  most 
aasilv  be  put  in  martingale  language  (see  Elliott, 
1982;  Kumar  and  Varaiya,  1986;  and  the  definitions 
below) : 


(0t,Vt)  in  (*)  «hould  be  a 


Assume  («t;  t6[0,TJ)  is  an  incrta?lnq  sequence  of 
information  fields,  i.s.  8s-lc0s;  *nF  •£(.!, T). 


A  random  saquencei  ((^)  is  said  to  be  a  Kartingale 
Diffarence  sequence  w.r.t.  fl t  *-tt  for  all  te[0,T], 


(1)  bt  1»  at_"***urab1*' 

(ii)  E(|it|)<-, 

(iil)  E(tt|8.)-0  «.*.»  for  all  s6[0,t-l). 


1 .  TIKE-REVERSION  OBJECTIVES 

We  want  to  obtain  a  time-reversed  version  of 

system  (1),  such  that  its  solution,  (?t>^t'0t)« 
in  some  sense  equivalent  to  (yt>xt<*tl*  To  make 
this  objective  explicit  it  needs  both  a 
specification  of  vhat  we  mean  by  a  time-reversion 
of  (1),  and  a  specification  of  the  desired  sense 
of  process  equivalence. 

By  a  reverse-time  system  we  mean  a  stochastic 
difference  equation  which  starts  at  time  T  and 
runs  in  negative  time  direction  on  the  interval 
(0,T).  We  require  from  a  time-reversion  of  system 
(1)  that  it  does  not  change  the  state  space  and 
that  the  solution  of  the  resulting  reverse-time 

system  represents  the  process  (pt,Rt,»t).  More 

speclficly,  (?t,9t,Jt)  must  be  the  solution  of  the 
following  system  of  stochastic  difference 
aquations,  all  t£(0,T-lJi 


Assuse  (»t;  tE(0,T])  is  a  decreasing  saquance  of 
information  fields,  l.e.  FsC9s.j;  any  s6(l,T). 

A  random  saquance  (Ct)  is  said  to  be  a 


It  iff  for  all  t£(0,T), 

(i)  tt  is  Sfmaasursbla, 

(ii)  E(|tt|)<-, 

(iil)  E((t|»s)-0  a.s.  »  for  all  s€(tel,T). 


Kav)~j  type  of  reverse-time 

system,  the  next  step  is  to  specify  the  types  of 
equivalence  of  solutions  of  systems  (1)  and  (4), 
in  which  we  are  interested.  For  stochastic 
processes  several  useful  types  of  equivalence  have 
been  defined  and  named  in  the  past.  Ws  restrict 
ourselves  to  the  two  most  Important  types  of 
equivalence  and  their  unambiguous  names  (Elliott, 
1982;  Jacod  and  Shlryaev,  1987); 

-  indistinguishable, 

'  -  equivalent  in  law. 

Definitions  are  given  below. 


1 _ Definition 


Two  processes  ( 1 1 )  an<l  <{t),  te(0,T],  ara  said  to 
b«  indistinguishable  if  they  ara  defined  on  the 
aaaa  probability  apaca  (a,»,P)  and 

P(  lt  -  tt  .  all  te(0,T)  >  -  1.  (5) 

4  Definition 

Two  procaaaaa  («t)  and  (lt|,  te(0,TJ,  ara  laid  to 
ba  aoulvalent  In  law,  if  thay  hava  tha  aaaa  -atata 
apaca,  K,  and  for  all  OitjrtjC _ <tkiT. 

P(U.  ...,r  )€dk)  -  P< <?  ...,I  J€dX)  ,  (6) 

C1  ck  C1  ck 

for  any  k  and  all  measurable  dXCE*. 

For  diacrata-tiaa  procasses  (5)  la  satiaflad  if 

and  only  if,  for  all  t£(0,T),  lt-Tt  *l«°»t  auraly. 
Our  objactlva  in  tha  sequel  la  to  obtain 
time-reversed  systems  of  typa  (4),  with  aolutiona 
that  ara  raapactivaly  Indistinguishable  and 
equivalent  in  law  w.r.t.  tha  solution  of  (l). 

4  INDISTINGUISHABLE  TIHE-REVERS ION 

In  this  section  wa  derive  a  typa  (4)  version  of 
ays tea  (1),  such  that  their  solutions, 

(?t>*t'Etl  *nd  (yt.*t.«tl.  ®r*  indistinguishable, 
and  illustrate  these  results  for  a  jump-linear 
example. 

The  first  step  of  our  derivation  consists  of  a 
substitution  of  (2)  and  (1)  in  (1),  to  arrive  at 
the  in  section  2  discussed  time-reversed  system, 

xt  *  »‘<»t+l.»f*t*l<vt>  <  (7. a) 

•t  "  b*(«t+l,vt),  (7.b) 

yt  -  c(st,xt,wt,utJ .  (7.c) 

Although  (7)  and  (4)  look  similar,  one  requirement 
is  not  met:  the  driving  noise  In  (7)  is  not  a 
reverse-time  Martingale  Difference  aequence  w.r.t. 
the  future  Information  field 

*t  s  •((ys,xs,«,,ws,vs,us);  se(t,T)|.  (8) 

Therefore  our  next  step  Is  to  Introduce  a 
particular  reverse-time  Martingale  Difference 
sequence,  (wt*,vt*|,  as  follows, 

<wt*»vt*)  "  («t«vt)  '  (“t<^t)  .  (9. a) 

with 

*t  »  E<vtl»t+ll«  (9.b) 

“t  *  E(wtl*t+l»vt»'  «U  t€(0,T-l) .  (9.c) 

and  (vT*,vT*)«o. 

Notice  that  the  definition  of  wt  differs 
significantly  from  the  reverse-time  predictable 
Process  ^(vtl^tilt*  1*  such  the  decomposition  in 
(9)  is  not  the  unique  canonical  decomposition  (see 
Jecod  and  Shiryaev,  1987).  The  introduction  of 
this  non-canonical  decomposition  is  a  crucial  step 
beceaeery  for  obtaining  the  time- ravers ion  of 
hybrid  state  system  (1). 

tn  the  sequel  we  verify  that  (wt*,vt*|  la  Indeed  a 
reverse-time  Martingale  Difference  sequence  w.r.t. 
»t,  and  thus  also  w.r.t.  9t*  a  u  «((*,*, v,*); 
•e(t,TJJ,  Moreover  we  show  that,  due  to  the 


duality  of  the  Markov  property,  (w,.,vt)  is 
conditionally  independent  of  *t+2  given 
(xt«-l>*t+l>  • 

5  Theorem 

Assume  (Wt,vt),  <wt,vt)  and  (wt*,vt*|  satisfy  (1) 
and  (9).  Then  (wt*,vt*)  is  a  reverae-tlme 
Martingale  difference  sequence  w.r.t.  9t*,  while 

wt  and  vt  satisfy: 

“t  “  E(wtl»t+l>»t'*t+l)«  (10. •) 

0t  -  *(vt|at<.1#xt+1|,  all  -v(0,T-l ) .  (10. b) 

Proof:  See  Blom  and  Bar-Shalom  (1989). 


Theorem 

S  implies  that  wt  and  vt  can  be 

written 

as 

o 

ft 

1 

f(t/®t+l'®t'xt+l) * 

(11 

■«) 

- 

g(t,*t+i.xt+l> • 

(11 

•b) 

Substitution  of  (9. a)  and  (11. a, b)  in  (l 

f.«,b,c) 

yields 

xt  " 

»{t,»t+l.®t.xt+l.w‘t)> 

(12. 

■“) 

®t  " 

h(t,St4.y.Xt+1,V*t)  , 

(12 

•  b) 

yt  - 

c(t,«tvl-»t-xt+l<xt-u‘t>ut>- 

(12. 

•  c) 

with. 

»(t,s,a, 

x,w*)  -  a*(s,a,x,w**f(t,s,a,x)), 

(13. 

a) 

b(t,s,x, 

,v*)  -  b*(S,v*+g(t,S,x)), 

(13. 

b) 

c(t,*,s, 

,x,z,w*.u)  -  c(a,z,w*+f (t,«,a,x) , 

u).  (13. 

c) 

The  above  result  is  summarized  by  the  following 
corollary. 

5 — cotellery 

Under  assumptions  A.l  and  A. 2,  the  solution 

(7t,SCt,»t)  of  the  reverse-time  system  (4)  is 
indistinguishable  from  the  solution  (yfxt’*t' 
system  (li  if 

( I )  (?T'*T'®t)  “  (yT'xT'*T^  *.*., 

(II)  I,  B  and  C  satisfy  (13.a.b,c), 

(III)  (0t,Pt)  "  (wt\vt*)  :  *11  t6(0,T-l} , 

with  w*t  and  v*t  satisfying  (9. a)  and  (10). 

To  illustrate  the  results  obtained  so  far,  let  us 
consider  the  particular  situation  of  a  linear 
system  with  first  order  Markovian  switching 
coefficients  and  observation  noise  independent  of 
the  system  driving  noise.  Both  a(S,a,x,w)  and 
c(a,x,w,u)  ere  then  linear  in  (x,w),  while  the 
first  is  a-invariant  and  the  second  is  w-invariant, 
by  which  system  (1)  simplifies  to, 
xt+l  “  *(#t4l)xt  +  B(*t+l)wt> 

•t+X  -  b(*t»vt)« 

Yt  "  c(*t>xt  ♦  K(«t>uf 
Then  from  Corollary  6  we  readily  find  the 
indistinguishable  time-reversed  system, 

xt  “  *-1(*t+l)  (xt»l  “  ®(*t+l>  (“t+w*t>)< 

'•t  ”  b*(®t+l,^t+v*t) « 
yt  -  C(*t)xt  +  H(St)ut, 


“here  <v*t,v*t)  is  the  r<v<ru-tlH  MD-sequanca  of 
Theorem  5,  .  vt 

end  f,  g  and  b*  are  according  to  (11)  and  (11.  b) . 
Tha  diffarence  aquation  for  xt  la  similar  to  tha 
ona  for  tha  linear  Cauaaian  axampla  in  aaction  2. 

But  dua  to  Ct,  it  may  avan  ba  nonlinaar  in.xt+1. 

At  tha  and  of  tha  naxt  aaction  va  will  show  that 
thara  ara  aoaa  furthar  aiapllficatlona  posaiMa 
for  thia  axaapla,  in  caaa  of  in  probability  law 
equivalence. 


5.  EQUIVALENT  IN  law  tike-reversion 


tl*t+i"*t'xt«-l 


(V|.)  - 


p«tl«t+l.*t.xt*i<“*v'tl*1 ' 


(It. a) 


vhara  vt  aatlafiea  (10. a). 

With  thia  our  remaining  atap  ia  to  characterize 
tha  danaity  at  tha  right-hand  aida  of  (14. a)  by 
applying  Bayea  formula. 


S  Proposition 

Undar  assumptions  A. 1  and  A . 4 .  tha  diatributlon  in 
(iv)  of  Thaoran  7  parmita  a  danaity  which  la 
characterized  by  (14. a)  and. 


In  thia  section  we  derive  conditions  under  which 
the  solutions  of  (1)  and  (4)  ara  equivalent  in  taw, 
and  discuss  these  results  for  a  jump-linear 
exampla.  So  far  our  line  of  reasoning  is  quite 
similar  to  tha  martingale  approach  of 
time-reversing  a  diffusion.  However,  things  ara 
quite  different  now  wa  require  equivalence  in  law 
only.  Tha  reason  ia  that  while  in  tha  diffusion 

situation  thia  raqulraa  that  dwt  and  dvt  ara 
equivalent  in  law,  no  similar  simple  results  hold 
in  the  discrete-time  situation.  Instead  of  this, 
ve  identify  tha  relation  between  conditional  laws 

of  wt  and  vt  by  a  Bayesian  eatimation  step.  Naxt 
we  characterize  f  and  tha  required  law  of  w*t. 


Under  assumption  A.l  tha  solution  (?t'*t>^t)  ot 
reverse-time  system  (4)  is  equivalent  in  law  w.r.t. 
the  solution  (yt,xt,«t)  of  system  (1)  if, 

(i)  PUJT.iT.ffxJedX)  -  PKyx.Xx.axlcdX); 

for  any  measurable  dXCX"xRnxJ4, 

(ii)  i  and  S  satisfy  (13. a, c) , 

(ili)  P(fft-4|fft+l-9.5tt+1-xJ  - 

-  P(9t-4(9t+1«9,xt+1«x), 

(iv)  P(0tedX|(St+1,»t+1,»t|-(x,«,a) )- 

-  P(wt*edX| (xt+l,at+l,»t)-(x,«,a) }, 
all  (x,9,4,t)eRnxHJx(Q,T-l)  and  measurable  dXCRP, 
with  w*t  and  f  satisfying  (9. a),  (10. a)  and  (11. a). 

Proof:  Sea  Blom  and  Bar-Shalom  (1989). 


Pvtl*t*l-»t.>tt*l 
.c(«,a.x)  p„ 

wt 


(.,a,a,x)  -  |vxa*T(s,a,x, .) 

(•)  Pv  (a*(a,n,x, .) (a) j, 
xtl*t 


(14. b) 


with  th«  gradient  *n<3  c  «ith«r  «  normalizing 
factor  or  zero  iff  p  (x|*,ii)«o. 


Moreover, 


P(»t-,'lSfH"*-*t+l_x)  ■  P(»r-n»fn“*)- 

•Pv  ,  (n | a , a )  p-1  (X|«).  (15) 


Proof:  See  Bloa  and  Bar-Shaloa  (1989) . 


Jump-linear  example 

For  a  linear  ayatem  with  first  order  Markovian 
switching  coefficients  we  arrived,  in  section  4, 
at  the  following  reversed-time  equation: 


xt  “  A-1  (9t+1)  (xt+l  "  B(*t<l)  [“t*"'tH' 
with  w\  the  reverse-time  MD  sequence  and 

wt-E{ wt I  * t+1 >  * t' xt+l * '  Because  a*  is  linear  in 
(x , w) ,  its  gradient  w.r.t.  x  is  w-invarlant,  by 
which  proposition  8  yields 
P. 

I  *t-+l  •  •  A  t  ♦  1 

A-l(9) (X-B(8)w) 


Wt'®t*l'#t'xt-H 

-  C}J9 , 4,  X)p 


(W|9,4,X) 

(«)P„ 


“tl't 


( A- 


4)  . 


In  spite  of  the  simplification  this  is  a  form 
which  is  in  general  quite  complex,  by  which  vt 
still  may  be  a  nonlinear  function  ot  xt4q. 
Obviously,  this  type  of  complexity  could  have  been 
expected,  as  it  is  well  known  that  a  discrete-time 
Bayesian  estimation  step  leads  to  nonlinear 
aquations,  unless  the  prior  densities  involved  are 
Caussian. 


Our  remalninq  problem  is  thm  chsracterization  of 
the  conditional  law  of  v*t.  As  this  is  actually  a 
discrete-time  nonlineer  filtering  problem,  it  can 
be  done  by  applyinq  Bayes  formula.  Me  do  this 
under  the  following  additional  assumptions: 

A.3-  The  a  priori  distribution  of  (x^,et) 
permits  a  density-mass  function  for  all  t6(0,TJ. 

AlxJL*  a*(9,4,x,w)  is  once  differentiable  in  xeRn 
for  all 

(»,4,w)€XjxRP. 

If  the  distributions  in  (iv)  of  Theorem  7  have 
density-mass  functions  then  It  can  easily  be 
verified  that  (iv)  Implies, 


6  CONCLUDING  REMARKS 

Ha  considered  the  problem  of  reversing  the  Markov 
solution  of  a  nonlinear  stochastic  difference 
equation  in  time.  The  nonlinearities  were  due  to 
nonlinear  coefficients  and  a  hybrid  state  space, 
i.e.  s  product  of  an  Euclidean  space  and  a 
discrete  set.  For  simplicity,  it  warn  assumed  that 
tha  process  in  tha  discrete  set  satisfies  the 
Markov  property.  Subsequently  ve  gave  a  precise 
description  of  our  time  reversion  objectives:  thm 
development  of  time  reversed  difference  equations, 
of  forms  similar  *o  the  original  aquation,  but 
driven  by  reversed-time  martingale 
difference  sequences,  such  that  thalr  solutions 
ara  respectively  indistinguishable  from  and  in 


probability  law  equivalent  to  the  eolutlon  of  the 
original  equation.  Following  thie  the  derivation 
Oi  the  indistinguishable  reverse-time  equation  wa* 
performed.  The  main  new  theoretical  reault  le  the 
introduction  and  evaluation  of  a  non-canonical 
(Jacod  and  Shiry*av,  1987)  reveree-time  martingale 
decomposition,  which  ie  appropriate  to  tae  hybrid 
etate  apace  situation.  In  contrast  with  this,  all 
previous  reverse-time  equations  are  based  on  a 
canonical  martingale  decomposition.  After  that,  it 
was  shown  how  the  in  probability  law  equivalent 
time  reversed  system  can  be  obtained  by 
Introducing  an  appropriate  5*y*sian  estimation 
etep.  As  expected,  this  Bayesian  estimation  step 
leads  to  closed  form  equations  whose 
dimensionality  often  complicates  further 
applications.  In  view  of  this,  in  Blom  and 
Bar-Sha)^m  (1989)  we  elaborate  the  Bayesian  step 
for  linear  systems  with  Markovian  switching 
coefficients  (jump-linear  systems),  and  apply  the 
the  results  to  smoothing  a  trajectory  with  sudden 
manoeuvers . 
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Abstract.  A  realistic  stochastic  control  problem  for  hybrid  systems 
with  Markovian  Jump  parameters  may  have  the  switching  parameters  In  both  the 
slate  and  measurement  equations,  fur  thermore.  both  the  system  stale  and  the 
lump  stales  may  not  be  perfectly  observed.  Currently  the  only  existing 
Implementable  controller  for  this  problem  Is  based  upon  a  heuristic  multiple 
model  partitioning  (MMP)  and  hypothesis  pruning.  In  this  paper  we  present  a 
stochastic  control  algorithm  for  stochastic  systems  with  Markovian  jump 
parameters.  The  control  algorithm  Is  derived  through  the  use  of  stochastic 
dynamic  programming  and  Is  designed  to  be  used  for  realistic  stochastic  control 
problems,  l.e„  with  noisy  slate  observations.  The  state  estimation  and  model 
Identification  Is  done  via  the  recently  developed  Interacting  Mutltple  Model 
algorithm.  Simulation  results  show  that  a  substantial  reduction  In  cost  can  be 
obtained  by  this  new  control  algorithm  over  the  (MMP1  scheme. 
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1.  INTRODUCTION 

An  Important  problem  of  engineering  concern 
Is  the  control  of  discrete-time  stochastic 
systems  with  parameters  that  may  switch  among  a 
finite  set  of  values.  In  this  paper  we  present 
the  development  of  a  controller  for  dtscrele-lime 
hybrid  lump-linear  Gaussian  systems.  Here  the 
state  and  measurement  equations  have  parameter 
matrices  which  are  functions  of  a  Markov 
switching  process.  The  j'lmo  states  are  not 
observed  and  only  the  stale  Is  observed  In  the 
presence  of  noise. 

Along  with  presenting  a  desirable  practical 
control  algorithm  we  also  point  out  an 
Interesting  theoretical  phenomenon.  We  show  that 
there  Is  a  natural  connection  between  the 
Interacting  multiple  model  (IMM)  slate  estimation 
algorithm  (Btl  and  the  control  of  jump-linear 
systems.  Thus  the  IMM  Is  the  stale  estimation 
algorithm  of  choice  for  use  In  these  types  of 
control  problems. 

Systems  which  pertain  to  the  Jump-linear 
modelling  methodology  are  found  In  many  areas. 
Systems  of  a  highly  nonlinear  nature  can  be 
approximated  by  a  set  of  linearized  models  (M3. 

VI.  V21.  A  failure  In  a  component  of  a  dynamical 
system  (or  subsequent  repalrl  can  be  represented 
by  a  sudden  change  In  the  systems  parameters  (82, 
Si,  VI).  Also  economic  problems,  which  can  be 
modelled  by  parameters  that  are  sub|ect  to  sudden 
changes  due  to  shortages  In  Important  materials 
(G2l.  And  as  Is  noted  In  (M6|  there  also  exist 
applications  to  the  design  of  control  systems  for 
large  flexible  structures  In  space. 

There  has  been  an  extensive  amount  of  work 
done  In  this  area  and  on  the  related  problem  of 
controlling  stochastic  dynamic  systems  with 
unknown,  time-invariant  parameters.  We  refer  the 
reader  to  the  (T3I  and  (G3j  for  a  list  of 
references  and  a  discussion  of  their  scope  and 
applications. 
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More  recently  In  (S21  a  feedforward/feedback 
controller  was  presented  for  the  continuous-time 
problem  with  a  completely  observed  system  state 
and  where  the  "modal  Indicator"  Is  measured  with 
a  high  quality  sensor.  In  (M61  the 
continuous-time  Jump-linear  problem  Is  considered 
where  the  system  state  and  "modal  processes^are 
perfectly  observed.  The  optimal  regulator  was 
obtained  and  notions  of  stochastic 
stablllzablllty  and  detectability  were  Introduced 
to  characterize  the  behavior  of  the  optimal 
system  on  long  lime  Intervals.  In  (H7)  the 
continuous-time  Jump-linear  problem  with  additive 
and  multiplicative  noises  and  noisy  measurements 
or  the  plant  stale  was  considered  with  the  plant 
mode  assumed  perfectly  observed. 

In  (Cl)  a  sufficient  stability  test  Is  given 
for  checking  the  asymptotic  behavior  of  the  error 
Introduced  by  the  averaging  of  hybrid  systems. 

In  (MSI  the  continuous-time  Jump-linear  problem 
with  non-Markovlan  ret  me  changes  was 
considered.  A  control  schearc  was  presented  for 
the  case  of  perfect  ovstrvatlons  of  the  system 
state  and  plant  regime. 

In  JC3J  a  discrete-time  Markovian  Jump 
optimal  control  problem  was  considered.  The 
controller  Is  for  the  case  of  perfect  system 
slate  observations  and  known  form  process.  They 
derive  necessary  and  sufficient  conditions  for 
the  existence  of  optimal  constant  control  laws 
which  stabilize  the  controlled  system  as  the  time 
horizon  becomes  Infinite.  Through  examples  they 
show  the  Interesting  result  that  stablllzablllty 
of  the  system  In  each  form  Is  neither  necessary 
nor  sufficient  for  the  existence  of  a  stable 
steady-stale  closed-loop  system. 

In  (YU  a  discrete-time  system  with  perfect 
state  and  mode  Information  was  considered.  A 
controller  was  presented  which  Is  stabilizing  In 
the  mean  square  exponential  sense. 

As  pointed  out  In  JG2J,  we  generally  cannot 
determine  the  optimal  Jump-linear  quadratic 
Gaussian  closed-loop  control  law  analytically 


even  ioi  a  iwo-sitp  problem.  In  order  lo  compute 
the  optimal  control  extensive  numerical  search 
methods  must  be  employed  and  thus  one  would  like 
to  find  slaipler  suboptlmal  control  schemes. 

Currently  the  only  existing  Implementable 
controller  f or  this  problem  (switching  parameters 
In  the  system  state  and  measurement  equations  and 
noisy  state  observations).  Is  the  one  discussed 
In  (T3)  and  Is  of  the  OLOF  class.  This  algorithm 
Is  based  upon  a  heuristic  multiple  model 
partitioning  (MMP)  and  hypothesis  pruning.  The 
MMP  approach,  being  simple  and  straightforward  to 
Implement.  Is  a  reasonable  choice  for  the  unknown 
parameter  problem  (Lll.  and  as  shown  in  (T31  It 
works  well  for  applications  Involving  switching 
parameters  In  the  state  measurement  equation 
only.  For  the  non-switching  parameter  problem 
the  operating  mode  Is  determined  lo  a  high 
probability  In  a  relatively  short  period  of  time 
and  the  MMP  approach  gives  the  linear  quadratic 
Gaussian  optimal  control. 

For  switching  parameter  problems  a  different 
situation  exists.  Here  because  of  the  switching 
the  operating  mode  may  not  be  determined  to  a 
high  probability.  The  proposed  approach  to 
deriving  a  suboptlmal  control  scheme  Is  to  start 
with  the  solution  to  the  optimal  control  problem 
via  the  use  of  stochastic  dynamic  programming. 

By  utilizing  dynamic  programming  and  making 
appropriate  suboptlmal  assumptions  the  use  of 
numerical  search  methods  has  been  avoided.  Ve 
thus  have  developed  a  multiple  model  control 
scheme  which  has  the  following  desirable 
properties:  (a)  It  gives  the  optimal  final 
control,  (bl  the  algorithm  utilizes  the  (MM  state 
estimation  scheme,  and  (c)  It  has  the  same 
property  as  the  MMP  approach  In  that  It  gives  the 
optimal  linear  quadratic  control  under  the 
assumption  of  a  perfectly  known  model  history 
sequence  (which  Is  however  an  unrealistic 
assumption  for  this  class  of  problems). 

For  comparison  purposes  we  Implement  the 
"switching  parameters  In  the  system  stale 
equation"  controller,  proposed  (but  not  tested) 

In  (T31.  Ve  show  via  example  that  a 
statistically  significant  reduction  In  cost  can 
be  achieved  through  the  use  of  our  controller 
which  also  belongs  to  the  OLOF  class. 

The  paper  Is  outlined  as  follows.  In  section 
2  the  problem  formulation  Is  given.  In  section  3 
an  Interesting  connection  between  the  IMM  state 
estimation  algorithm  and  the  control  of  multiple 
model  systems  Is  shown  to  exist.  In  section  1  we 
obtain  the  control  algorithm.  A  new  "full-tree" 
control  algorithm  Is  derived  which  utilizes  all 
possible  future  parameter  history  sequences.  In 
section  5  we  use  simulations  to  compare  the  MMP 
control  algorithm  with  the  full-tree  controller. 

2.  PROBLEM  FORMULATION 

The  problem  to  be  solved.  Is  discussed  next. 

Ve  look  the  pragmatic  approach  of  starting  with 
the  available  mathematical  and  statistical  tools 
found  to  yield  success  In  solving  similar 
problems  of  this  type  In  the  past  (l.e.,  use  Is 
made  of  the  stochastic  dynamic  programming  method 
and  the  total  probability  theorem,  etc.).  As  we 
shall  see,  not  only  does  this  practical 
engineering  approach  yield  an  Improved  multiple 
model  control  algorithm,  but  It  also  leads  to  the 
interesting  theoretical  observation  of  a  direct 
connection  between  the  IMM  state  estimation 
algorithm  and  jump-linear  control. 

It  Is  desired  to  find  a  sequence  of  causal 
control  values  to  minimize  the  cost  functional 


J  *  e(c(0))-e(x(N|0(N)x(N).£  |x(k)-Q(k)x(k) 

•  u(k|'R(k)u(k)|j  (2.1) 

where  Q(k|tO  for  each  k»0.1....N  and  and  It  Is 

sufficient  that  R(k)>0  for  each  k*0,l . N-t. 

The  discrete-time  system  state  and 
measurement  modeling  equations  are 

xfk)  -  F(M(k)|x(k-l)  •  G[M(k)]u(k-l) 

•  v(k-l,M(k )  |  (2.2a) 

zlk)  -  H(M(kJ|x(k>  ♦  w(k.Mlk))  k-0,1.2.  ..  ( 2.2b) 

where  xfk)  Is  an  nxl  system  state  v  ctor,  u(k)  Is 
an  pxl  control  Input,  and  z(k)  Is  an  mxl  system 
state  observation  vector.  The  argument  M(k| 
denotes  the  model  "at  time  k"  -  In  effect  during 
the  sampling  period  ending  at  k.  The  process  and 
measurement  noise  sequences,  v(k-l,M(k|]  and 
w(k.M(k)l,  are  while  and  mutually  uncorrelaled. 

The  model  at  time  k  Is  assumed  to  be  among  a 
finite  set  of  r  models 


Mlk)  t  (1.2 . r) 

(2.3) 

for  example 

FlM(k)-j)  -  Fj 

(2  1) 

v(k-!.M(k|*J|  '  A’ttr.V,! 

(2.5) 

w|K.M(k)-Jl  -  7/(k  V  1 

(2.6) 

l.e..  the  structure  of  l hi  system  and/or  the 
statistics  of  the  noises  might  be  different  from 
model  to  model. 

The  model  switching  process  to  be  considered 
here  Is  of  the  Markov  type  The  process  Is 
specified  by  a  transition  matrix  with  elements 

P.j-  Lel 

1*  £  (z(O).z(l) . zlkl.ulOl.ulll . ulk-l ) I  (2  7) 

denote  the  Information  available  lo  the 
controller  at  time  k  (l.e.  the  control  Is 
causal). 

3.  THE  LAST  STAGE  CONTROL  AND  THE  CONNECTION 
WITH  THE  IMM  ESTIMATOR 

An  Integral  part  of  any  control  algorithm  for 
this  class  of  problems  is  the  system  state 
estimator.  In  this  section  we  show  that  there 
exists  an  Interesting  connection  between  the 
control  of  multiple  model  stochastic  systems  and 
the  IMM  system  state  estimator  (Bll.  To  this  end 
we  start  by  solving  for  the  lime  N-l  optimal 
control.  The  optimal  control  at  lime  N-l.  Is  the 
vaiuz  of  u(N-l)  which  minimizes 

JlN-1)  -  E{x(N-l)  Q(N-l)x(N-l).u(N-l)'R(N-l)u(N-l) 
♦x(N|'Q(N)x(N)|l',",| 

-  £  E{x(N-l)’Q(N-l)x(N-l)*u(N-irR(N-l)u(N-ll 
yi 

♦x(N)'Q(N)x(N||l"'1.M(N).j) 

•  P(M(N|-)||‘''1)  (31) 


«  p(h(ni*jiin',j  (3  2) 

and  use  the  state  equation  (2.2al  and  (2.-1) 

(2.5)  In  13.1)  to  get 

J(N-l)  -  I  c{x(N-l)'|o(N-l  I*Fj'Q(N)fJx(N-1  | 

♦2u(H-l)'Cj'Q(N)rjxlN-l)»u(N-l)'[R(N-l)*CJ'Q(N)Cj] 

•u(N-I)|(n',.M(N).j}mj|N|N-1) 

♦£  trlQlNlV.IUjINlN-l)  (3.3) 

Now  taking  the  partial  of  (3.3)  w.r.t.  u(N-l)  and 
setting  It  to  zero  yields 

u'(N-i)  -  -Ir(n~i)»£  g’qihig.m .( N|N— I ) ]  ‘ 

1  )=i  1  11  1 

•  £  G'.Q(N)F.E{x(N-l|lr\hlN)-j} 
yi  1  ’  1 

•P^NIN-I)  (3.1) 

Notice  that 

e{x(N-1)||"'1.M(Nm)  -  £  E{x(N-l)|lN'1.M(N)-i, 

I  .51 

M(N-1)-|}  P|M(N-l)-l|M(N).l.rMl  (3.5) 

where,  since  M(N)*J  In  the  first  conditioning  Is 
Irrelevant,  the  expectation  Inside  the  summation 
Is 

e{x(N-1)|iN1.M(N)"jJ  -  £  XjlN-llN-DPy.lN-llN-l) 

£  x0i(N-l|N-l)  (3.S) 

which  Is  the  IMH  mixed  Initial  estimate  (81). 

Thus  using  (3.G)  In  (3.1)  we  gel 

u’(N-l)  •  *[R(N-l)*£GjQ(N)p^(NtN-l))"1 

•  £G'.Q(NlF.x0i(N-l|N-l)p  (NIN-I)  (3.7) 

Pi  1  1  1 

1.  THE  CONTROL  ALGORITHM 

We  will  derive  a  full-tree  control  algorithm 
(FT)  which  computes  control  values  by  taking  Into 
account  all  possible  future  model  histories.  As 
will  be  seen  by  our  example  this  method  offers 
Improved  performance  over  the  existing  scheme 
I T3 ). 

The  l-th  future  history  of  modtls  Is 
denoted  as 

-  (M(k)-I, . M(N)-IH)  l-l.....r"'l*‘  MU 

where  I,  Is  the  model  at  time  I  from 
history  I  and 

1  s  I,  i  r  l-k . N  (1  21 


J’(k.ll|  5  mljt  c{x(k)'Q(k)x(k)«u(k)*R(k)u(k) 

*  J-(k.l.l‘Ml|l‘}  (1.31 

where  J*(k.l‘|  is  the  optimal  cost-to-go  from 
time  k  to  the  end.  Now  applying  the  total 
probability  theorem  to  (1.3)  yields 

J"(k,lk)  •  mtn  I  (e( x(k )‘Q(k ) x(k)  ♦  u(k)'R(k)u(k) 
x(0  |,l  1  1 

♦  J'(k*l.l‘*,)|M‘*,-"J.l‘}p|M*'l-NJ|l‘))  (1.1) 

The  control  that  minimizes  an  approximation 
to  Hi)  Is  derived  In  the  Appendix,  and  Is  given 
as 

rM'k#2  -i 

ufT(k)  ■  -  [fifkl  ♦  I  G^pWllG^ujNlkMl] 

•  I  G,'  P‘(k*I)f.  x“(k|k)p  (N|k*l|  (1.5) 

i,i  ***  '*> 

and  again  we  see  the  natural  way  the  IMh  mixed 
Initial  estimates  show  up. 

Note  that  the  control  parameters  p'lk) 
(model-history-condltloned  optimal  cost  matrices) 
are  computable  off-line. 

5.  .SIMULATION  RESULTS 

The  FT  controller  developed  In  Sec.  1  Is  used 
to  control  the  stale  trajectory  of  the  system. 

The  performance  or  this  algorithm,  as  determined 
by  (2.1),  Is  compared  to  the  cost  obtainable  by 
using  the  MMP  controller  discussed  In  (T3J.  In 
order  to  obtain  a  meaningful  comparison  we  use 
the  rigorous  statistical  analysis  technique 
presented  In  [85.  W3J. 

The  control  of  a  double  Integrator  system 
with  process  and  measurement  noises  Is  considered 
with  a  gain  failure.  The  two  possible  models  are 
given  by  the  following  system  equation 


nr 

1  x‘(k)  . 

0 

u(k) 

L  o  ,  . 

J 

.  b*  . 

J 

J  v(k )  1*1.2  (5.1) 

with  measurement  equation 

z(k)  •  (I  0)  x‘(k)  ♦  w(k)  (5.2) 

The  models  differ  In  the  control  gain  parameter 
b‘.  The  process  and  measurement  noises  are 
mutually  uncorrelated  with  zero  mean  and 


variances  given  by 

E(v(k)  v(J))  -  0.16  «kj 

(5.3) 

and 

E(w(k)  w| J)  1  *  6tJ 

(5  1) 

The  control  gain  parameters  were  chosen  to  be 
b‘*2  and  b’-O.S. 

The  Kartkov  transition  matrix  was  selected  to 
be 


(5.51 


L  0.1  0.9  J 

For  this  example  N-7,  and  ihe  cost  parameters 
R(k)  and  Q(k),  (see  (2.1)1,  were  selected  as 

R(k)  •  5.0  k-1,2 . N-l  (5.6) 


the  FT  controller  performs  belter  than  the  MMP 
controller  (or  this  problem.  The  estimated 
Improvement  (decrease  In  cost)  of  70X  Is 
statistically  significant. 


and 
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_  _j 

J°d°  8:8. 

(5.7) 


where  the  last  matrix.  Q(7),  reflects  our  desire 
to  drive  x,(7)  vigorously  to  zero.  Also  note 
that  for  this  example  T-1.0. 

The  real  system  was  Initialized  with 
x(0)*(30.0.  0.01'  and  a  random  selection  was  done 
for  choosing  the  Initial  model  with 
P(M(0)*I]*0.5,  l»l,2.  The  Kalman  filters  each 
received  an  Initial  state  covariance  of 


P(0,0)  ■  [  1.0  2.0  ] 

and  the  Initial  stale  estimate  was  selected  as 


r  X|(oioi  “1 

r  z(o1  i 

L  XjColQ)  J 

L  z(0)  -  Z(-l]  J 

where  z(*tl  ■  30.0  ♦  w(-|)  and  z(0)  «  30.0  » 
w(0). 

Statistical  tests  were  made  on  the  results  of 
50  Monte  Carlo  runs.  Sample  means  and  variances 
-  of  the  Monte  Carlo  costs  C(  defined  In  (2.11 
were  computed  for  the  FT,  MMP,  and  "known 
mcdel-htstory"  (l.e.  optimum  linear  quadratic) 
controllers. 

Table  I  contains  the  results.  The  FT 
algorithm  shows  a  clear  reduction  In  cost  as 
compared  with  the  MMP  scheme.  However  In  order 
to  provide  a  rigorous  argument  that  the  actual 
performance  Is  ordered  as  Table  I  Indicates  we 
apply  the  statistical  test  presented  In  (BS,  V31. 

Table  II  contains  the  results.  The  sample 

standard  deviation  o-  of  the  mean  of 
the  cost  differences,  c“Mf-cP,  are  shown. 

The  hypothesis  that  the  FT  controller  Is  better 
than  the  MMP  scheme  can  be  accepted  only  If  the 
probability  of  error  a  Is  less  than,  say,  1 
percent.  Then  the  threshold  against  which  we 
compare  the  test  statistic  A/Oj  Is 
P"2.33.  This  test  statistic  has  to  exceed  the 
threshold  In  order  to  accept  the  hypothesis. 


TABLE  I 


TABLE  II 

STATISTICAL  TEST  FOR  ALGORITHM  COMPARISONS 


Test 

Estimated 

Statistic 

Improvement 

A 

°i  l/oi 

r 

FT-MMP 

13,956 

3,316  9.1 

70 

6.  CONCLUSION 

The  development  of  a  new  control  algorithm 
for  discrete-time  hybrid  stochastic  systems  with 
Markovian  Jump  parameters  has  been  presented. 

This  contoller  was  derived  through  the  use  of 
stochastic  dynamic  programming  and  by  taking  Into 
account  alt  possible  future  "histories  of 
models".  This  scheme  uses  the  1MM  state 
estimation  algorithm  Ve  show  that  there  Is  an 
Interesting  connection  between  the  IMM  state 
estimator  and  control  of  Jump-linear  hybrid 
systems.  This  new  controller  Is  of  the  DL0F 
class  and  has  off-line  computable  control  gain 
parameters. 

From  the  example  It  Is  seen  that  this  scheme 
can  achieve  a  statistically  significant  reduction 
In  cost  when  compared  to  the  multiple  model 
partitioning  approach. 


APPENDIX 


I.  Derivation  of  (1.5) 

Note  that  given  the  future  history  of 
models  m‘,lmj,  the  optimal  cost-lo-go 
J‘(k»l.lk,ll  Is  easily  computed  and  Is 
denoted. 

Jj(k-*!.l‘*,|  £  E{x(kMI'Pl(kM)x(k«l)|lk“,Mk*u'J) 

♦  oc’(  k*  I )  (A.l) 

where  the  notation  from  (B9)  Is  used  for  P(k»l) 
and  a(k*l). 

Since  the  expectation  In  (9.9)  Is  conditioned 
on  Mk,u<J,  we  obtain  our  of  approximation 
by  replacing  J*(k»l,lkM)  inside  the 
expectation  with  (A.l),  and  (9.91  becomes 


SAMPLE  AVERAGE  COSTS  AN0  STAN0AR0  0EVIATI0NS 


Known 

Mo  del -His  lory 

FT 

MMP 

Sample  Mean 

2,697 

6,063 

19,519 

Sample  Standard 
Deviation 

8,096 

3.96E5 

1.12E7 

rM-k«J 

J‘(k.lk)  -  mlji  X  |c{x(k)'Q(k)x(k)  ♦  u(k)'R(k)u(k) 

♦  c[x(k»l)*  Pl(k*l)x(k*I)|lk*l,Mk'u,JJ 

♦  o'(k.|)|Mk'lx,,lk)li||N|k.l|)  (A.2) 


where 


U,(N|k«l|  feP(fl‘*u,J|l‘|  (A.3) 

Now  use  (2.2a)  and  apply  the  smoothing 
property  of  expectation  to  (A. 2)  to  get 

fX-k*! 

J'(K.I*)  «  min  I  [c{x(k)'Q(k)x(k)  »  u(k)'R(k)u(k) 

i.i  '  ' 

♦  [F^^lkl'C^ulkl  ♦  v[k-Uk.,|]  P’(k*l  )[•] 

♦  ot,Ck*ll|Mfc*1-WJ.lk}n1CN|k*l|J  (A.4) 

Take  the  partial  w.r.t.  u(k)  of  (A. 41  and  set  to 
zero  to  solve  for 

rK-k«2 

u*(k|  —  [fi(k)*  £  Cik.IP,(k»l|CltMUI(N|k*l||'1 

*  I  Gl'k.lP'(kM)FlkME{x(k)|Mk4l-MJ.|‘}u|(N|k.l|  (A.5) 

We  still  need  to  evaluate  the  expectation  In 
(A.5).  This  Is  done  as  follows.  Note  that 

x(k)  Is  Independent  of  MCI).  l-k»2 . N  If 

M(k*ll  Is  known,  thus 

e{x(k)|n‘*1>J.l‘}  -  E{x(k)jf1(k*l)  -  lk.,.l‘)  (A. 6) 

But  (A. 6)  Is  x“(k|k).  the  IMM  mixed 
Initial  estimate  (see  (3.6)).  thus  using  (A.6)  In 
(A.5),  we  get  (4.5). 
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ABSTRACT 


Piecewise  Deterministic  (PD)  Markov  processes  form  a 
remarkable  class  of  hybrid  state  processes  because, 
in  contrast  to  most  other  hybrid  state  processes, 
they  include  a  jump  reflecting  boundary  and  exclude 
diffusion.  As  such,  they  cover  a  wide  variety  of 
‘-'pulsively  or  singularly  controlled  non-diffusion 
processes.  Because  PD  processes  are  defined  in  a 
pathwise  way,  they  provide  a  framework  to  study  the 
control  of  non-diffusion  processes  along  the  same 
lines  as  that  of  diffusions.  An  important 
generalization  is  to  include  diffusion  in  PD 
processes,  but,  as  pointed  out  by  Davis,  combining 
diffusion  with  a  jump  reflecting  boundary  seems  not 
possible  within  the  present  definition  of  PD 
processes.  This  paper  presents  PD  processes  as 
pathwise  unique  solutions  of  an  ltd  stochastic 
differential  equation  (SDE) ,  driven  by  a  Poisson 
random  measure.  Since  such  an  SDE  permits  the 

inclusion  of  diffusion,  this  approach  leads  to  a  - 

large  variety  of  piecewise  diffusion  Markov 
processes,  represented  by  pathwise  unique  SDE 
solutions. 

JL _ INTRODUCTION 

Because  many  of  the  stochastic  processes  that  we 
meet  in  nature  have  a  state  space  that  is  a  product 
of  a  continuous  space  and  a  discrete  set,  we  often 
need  pathwise  models  on  such  a  hybrid  state  space. 

As  a  result,  several  classes  of  hybrid  state  space 
models  have  been  developed,  such  as  systems  with 
Markovian  switching  coefficients,  doubly  stochastic 
counting  processes  and  Markov  decision  drift 
processes.  These  models  are  used  in  quite  different 
fields  of  applications,  by  which  their  studies  have 
often  evolved  separately.  One  reason  to  study  hybrid 
state  space  processes  within  a  common  framework  is 
that  their  martingale  parts  are  in  general 
discontinuous.  This  property  has  attracted  a  lot  of 
attention,  and  is  by  now  very  well  documented 
(Jacod,  1979;  Cinlar  et  al.,  1980;  Bremaud,  1981; 
Elliott,  1982;  Bensoussan  and  Lions,  1984;  Ethier 
and  Kurtz,  1986;  Jacod  and  Shiryaev,  1987).  It  is 
quite  clear  from  these  results  that,  to  study  hybrid 
state  Markov  processes  along  the  same  lines  as 
diffusions,  we  need  both  pathwise  representations 
and  strong  Markov  (martingale)  characterizations  of 
those  processes.  Unfortunately,  for  hybrid  state 
Markov  processes  there  is  presently  a  lacuna  of 
pathwise  representations  with  strong  Markov 
characterizations.  This  lacuna  is  apparent  if  we 
depict  the  main  classes  of  hybrid  state  Markov 
processes  in  the  form  of  a  Venn-diagram  (fig.  1) . 


There  exist  pathwise  representations  with  strong 
Markov  characterizations  of  counting  processes  with 
diffusion  intensity  (Snyder,  1975;  Marcus,  1978),  of 
diffusions  with  Markovian  switching  coefficients 
(Wonham,  1970;  Brockett  and  Blankenship,  1977)  and 
of  Piecewise  Deterministic  (PD)  Markov  processes 
(Davis,  1984) .  For  many  other  Markov  processes  in 
figure  1,  there  exist  only  strong  Markov 
characterizations  (Kingman,  1975;  Anulova,  1979, 
1982;  Bensoussan  and  Lions,  1984;  Belbas  and 
Lenhart,  1986) .  Actually,  PD  Markov  processes  seem 
the  most  interesting  of  all  processes  in  figure  1, 
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as  they  provide  pathwise  representations  with  a 
strong  Markov  characterization  of  all  major  non¬ 
diffusion  Markov  processes.  As  such,  PD  Markov 
processes  provide  *  irumcwork  to  study  Markov 
decision  drift  processes  (Hordijk  and  Van  der  Duyn 
Schouten,  1983;  Yushkevich,  1983)  along  the  same 
line  as  diffusions  (Vermes,  1985) .  With  this,  an 
interesting  generalization  is  to  extend  the  spectrum 
of  hybrid  state  Markov  processes  by  including 
diffusion  into  PD  Markov  processes.  As  the  present 
definition  of  PD  processes  does  not  seem  to  have  an 
opening  left  for  that  inclusion  (Davis,  1984),  we 
need  a  different  approach. 


Fig.  1.  Main  classes  of  hybrid  state  Markov 


processes. 

The  approacfiHthat  overcomes  this  difficulty, 
presented  in  the  sequel.  Is  to  assume  a  stochastic 
differential  equation  (SDE)  in  a  hybrid  space  and  to 
construct  a  large  class  of  piecewise  diffusion 
Markov  processes  from  it.  With. respect  to  the  state 
space  we  restrict  our  attention  to  a  hybrid  subset 
of  a  Euclidean  space.  Then  the  most  general  SDE  is 
of  Its  type,  driven  by  Brownian  motion,  w,  and  a 
Poisson  random  measure,  p  on  (0,“)xU, 

d*t  “  «(et>dt  +  0(!t)dwt  +  i  *(tt-'uJ  p(dt,du)  . 

The  path  of  a  solution  of  this  SDE  is  right 
continuous  and  has  left  hand  limits:  tt_  -  JJg  *t-A- 

If  p  generates  a  multivariate  point  (t,ut),  then  the 
path  of  t  has  a  discontinuity: 

*t  “  *t-  +  *»t-<ut>  • 

In  the  sequel  we  shall  focus  on  pathwise  unique 
solutions.  The  classical  result  for  the  existence  of 
such  solutions  requires  that  *  is  sufficiently 
continuous  (Gihman  and  Skorohod,  1972),  which 
restricts  the  SDE  essentially  to  systems  with 
Markovian  switching  coefficients.  However,  there  are 
some  non-classical  pathvise  uniqueness  results  that 
allow  a  discontinuous  *  (Lepeltier  and  Marchal, 

1976;  Jacod  and  Protter,  1982;  Veretennikov,  1988). 
Taking  these  results  as  a  starting  point,  we 
introduce  and  evaluate  a  particular  structure  for  * 
in  section  2.  This  structure  poses  hardly  any 
restrictions  on  the  possible  solution  of  the  SDE, 
while  it  enables  a  separate  evaluation  of  an 
unbounded  jump  intensity  and  a  hybrid  state  space 
situation.  In  view  of  this  separation,  we  first 
consider,  in  sections  3  and  4,  the  modelling  of  a 
jump  reflecting  boundary  in  through  an  unbounded 
jump  intensity,  and  after  that,  in  section  5,  we 
consider  the  hybrid  state  situation. 

Assume  an  open  subset  0  of  Rn  with  jump  reflecting 
boundary  #0,  which  means  that  (tt)  undergoes  an 


instantaneous  jump  into  the  Interior  of  0  if  (it;) 
tries  to  cross  or  to  travel  through  30.  To  model 
this  with  the  above  SDE,  the  Poisson  random  measure 
p  should  Instantaneously  generate  a  point  when  (Et) 
enters  30.  However,  this  is  not  possible  as  a 
Poisson  random  measure  generates  almost  surely  no 
point  at  an  entrance  time.  To  overcome  this  problem, 
we  briefly  discuss  the  following  three  approaches: 

1.  Replace  p  by  a  random  measure,  with  almost 
surely  one  point  at  an  arbitrary  time. 

2.  Assume  a  t  such  that  p  generates  an  active  point 
during  an  infinitesimal  small  time  interval 
after  entering  30. 

3.  Assume  a  t  such  that  p  generates  an  active  point 
during  an  infinitesimal  small  time  interval  just 
before  entering  30. 

Approach  1  adequately  solves  the  instantaneous  jump 
problem  but  creates  many  new  problems,  because  if  p 
is  not  a  Poisson  random  measure,  then  the  SDE  can 
not  be  analysed  within  the  powerful  It6  framework. 
Approach  2  Is  the  well  known  approach  of  randomized 
stopping  (Bensoussan  and  Lions,  1984).  As  this 
approach  allows  (it)  to  cross  or  to  travel  through 
■’O  thi  resulting  pi’  jc.ss  it  at  bet.  a  modnicataon 
of  a  PD  Markov  process.  Approach  3  is  the  desired 
solution.  However,  the  problem  with  approach  3  is 
that  it  is  in  general  not  known  how  to  carry  it  out. 
A  constructive  answer  to  this  will  be  given  in  the 
sequel.  It  is  clear  that  approach  3  needs  a  kind  of 
prediction  of  the  time  that  (tt)  might,  otherwise, 
enter  30.  Actually,  PD  Markov  processes  are 
presently  the  only  processes  for  which  this 
prediction  problem  is  solved  (Davis,  1984).  As  such, 
we  first  formulate  that  solution  in  an  SDE  set  up  in 
section  3.  Next,  in  section  4,  we  present  a  solution 
of  the  prediction  problem  for  the  situation  with 
diffusion . 

Finally,  in  section  5,  we  explicitly  consider  the 
hybrid  state  space  situation.  The  most  interesting 
effect  of  the  hybrid  state  space  assumption  is  that 
it  leads  to  a  particular  type  of  jumps:  jumps  in  the 
continuous  state  component  of  (Et)  that  anticipate  a 
simultaneous  transition  of  the  discrete  component  of 
(tt).  This  type  of  jumps  have  been  introduced  by 
Gnedenko  and  Kovalenko  (1968)  for  piecewise  linear 
processes  and  by  Sworder  (1972)  for  systems  with 
Markovian  switching  coefficients.  For  short  we  refer 
to  these  anticipating  simultaneous  jumps  as  hybrid 
lumps ■  The  SDE  framework  of  this  paper  provides  an 
elegant  way  of  representing  the  hybrid  jumps  of  PD 
Markov  processes  and  their  piecewise  diffusion 
generalizations. 

Some  other  interesting  generalizations  of  pd  Markov 
processes,  not  considered  in  the  sequel,  are  the 
inclusion  of  continously  reflecting  or  sticky 
boundaries.  The  inclusion  of  a  continuously 
reflecting  boundary,  while  preserving  pathwise 
uniqueness,  seems  possible  if  that  boundary  is 
smooth  enough  (Chaleyat-Maurel  et  al.,  1980;  Menaldi 
and  Robin,  1985;  Frankowska,  1985;  Saisho,  1987) . 

The  inclusion  of  a  sticky  boundary  without  loosing 
pathwise  uniqueness  seems  difficult  if  not 
impossible,  but  strong  Markov  characterizations  are 
possible  (Kingman,  1975;  Anulova,  1979,  1982). 


R+  -  (0,®)  and  R~  -  (-®,0), 

R+  «  R  +(0)  and  R_  «  R~+(0). 

I  -  (..,-2, -1,0, 1,2,..), 

N  -  (1,2,3, ..) . 

1  -  Col((2 - 5n)  if  *  „  Col( 1 1(  , 

l*l2  =  "ii2  »  if  «  is  a  matrix 

1*1  “  J  »i  ,  if  *  is  s  vector 


30 
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•  i-th  component  of  process  Et* 
boundary  of  the  closure  of  set  0. 

:  integer  part  of  x. 

:  x (True) -1  and  X(False)-0. 

:  right  continuous  with  left  hand  limits 
:  the  set  of  all  real-valued  functions  that 
are  k  times  continuously  differentiable  on 
0.  The  superscript  is  deleted  if  k-o.  If  k 
is  followed  by  b,  then  f  and  its  first  k 
derivatives  are  bounded  on  0. 

:  domain  of  operator  4. 


We  assume  a  stochastic  basis  (0,»,r,P),  endowed  with 
an  m-dimensional  standard  Wiener  process,  (w.),  and 
a  Poisson  random  measure,  p(dt,du)  on  Rj.xRd+*  (Jacod 
and  Shiryaev,  1987,  p.  70),  with  intensity  measure 
dtxm(du) ,  and  consider  the  following  stochastic 
differential  equation  (SDE)  in  R+xRr‘, 
dt t  -  a(tt>dt  +  fl(it)dwt  +  R-'Rd  *at-<u>  q(dt,du)  + 

+  R+xRd  *(tt->u)  P(dt,du)  ,  (1) 

where  q  is  the  martingale  measure  of  p,  Eq  is  an 
Sg-measurabl  -  random  variable,  while  a,  B  and  *  are 
measurable  mappings  r'r  appropriate  dimensions. 

The  classical  reference  for  equation  (1)  is  Gihman 
and  Skorohod  (1972).  Significant  exLtr.'ions  of  their 
results  have  been  obtained  by  Lepeltier  ana  Marchal 
(1976)  in  their  study  of  the  relation  between  an 
integro-dif ferential  operator  and  an  SDE  of  type 
(1).  Their  particular  SDE  can  easily  be  obtained 
from  (1),  by  introducing  homeomorphism  mappings  of 
R~XRd  jnto  (u£Rd+1;0<|u| £1)  and  of  R+xRd  into 
(uti»d'1 it  |  u|  <•) ,  anu  subsequently  transforming  m 
and  correspondingly.  Consequently,  the  results  of 
Lepeltier  and  Marchal  can  immediately  be  used  in  the 
present  study  of  (1),  while  allowing  the  intensity 
of  the  active  points  in  R+  to  be  unbounded  outside 
some  known  Borel  set  0'CRn. 


A ■ 1  There  is  a  constant  K  such  that,  for  all  E6Rn, 
I  a (C )  | 2  +  IS  (t ) I2  +R_/Rd  l*(E,u)  |2m(du)  £  K(1+|EI2). 

A. 2  For  all  kEN  there  exists  a  constant  La  such 
that,  for  all  E  and  y  in  the  ball  Bk=(uERn; |u|2£k) , 

|  a  (E  )  -a  (y)  |  2  +  |B(E)-fi(y)l2  + 

+  R-£Rd  l*(E.u)-*(y,u)  |2m(du)  £  LjJ  E  -y  |  2  . 

A'  .  3  0'  is  a  known  Borel  subset  of  Rn, 

R+xRd  * ^  *(E,u)*0  )m(du)  is  uniformly  bounded  on  O', 

and  [;+*•(!, u)]  6  O',  for  all  EeRn,  uERd+1. 

A' ■ 4  For  all  kEN  there  exists  a  constant  K^,  such 
that,  with  Bj.  the  ball  of  A.  2 : 
a.  for  all  E6Bkno', 

R+xRd  l*(E'u)  1  m(du)  S  Mk- 
fc.  for  all  EEBkn(Rn-0') , 

^  £d  l*U.u)  I  m(du)  £  Mk, 

given  that,  for  all  uER+xRd, 

*(E,u)  -  *(E,u+Col(l,0,..,0)) . 

A ' ■ 5  For  all  rEN  there  is  a  constant  NT,  such  that 

E(  l  R+'R d  * (  *(Es-,u)*0  )  p(ds,du))  £  Nt . 


Given  m(du) “du^xx (du)  and  assumptions  A. 1 ■  A. 2 . 

A' . 3 ■  A' ■ 4 ■  A'  ■ 5  are  satisfied.  Then  equation  (1) 
has  for  any  Eoeo'  a  pathwise  unique  solution,  (Ep). 
Moreover  (Et>  is  then  a  right  continuous  Markov 
process. 

Remark:  Proposition  2.1  is  a  version  of  Theorem  III4 
of  Lepeltier  and  Marchal  (1976),  in  the  sense  that 
they  considered  the  situation  of  O'-  Rn. 
Nevertheless,  for  the  proof  we  can  almost  follow 
Lepeltier  and  Marchal.  Another  recent  extension  of 
Theorem  III4  of  Lepeltier  and  Marchal  is  to  the 
situation  of  a  non-Lipschitzian  a  in  turn  of  a 
sufficient  non-degeneracy  assumption  on  B 
(Veretennikov,  1988)  . 

E£22l: 

If  (l)'s  fourth  right  hand  term  vanishes,  then  it  is 
well  known  that  A ■ 1  and  A ■ 2  are  sufficient 
conditions  (Gihman  and  Skorohod,  1972)  .  As  such,  we 
have  to  show  that  (l)'s  fourth  right  hand  term  does 
not  change  that  situation,  under  A' . 3 .  A' ■ 4  and 
M.5- 

Due  to  A' . 3  and  the  definition  of  It6  integration  a 
solution  of  (1)  is  CADLAG.  Due  to  A' . 5.  the 
discontinuities  in  (Et).  that  are  caused  by  (l)'s 
fourth  right  hand  term,  are  countable.  Therefore  we 
can  associate  with  each  discontinuity  a  time,  T4, 


such  that 


arvT~a  multi-variate  point,  u  , 

Td 

0<T1<T2<. .<Tj<. .  and  |ia  Td  «  ®. 
and  (Et)  being  CADLAG, 


Due  to  the  latter 


the  latter  sum  is  finite  (a.s.)  for  all  t£Rj.,  due  to 
A' ■ 4  and  A' ■ 5 ■  With  this  result  it  is  sufficient  to 
shew  that  (1)  has  a  pathwise  unique  solution  on  an 
arbitrary  finite  time-interval  (0,T).  For  the 
existence  of  a  solution,  see  the  proof  of  Th.  Ill* 
of  Lepeltier  and  Marchal  (1976;  pp.  82-85).  We 


already  know  that  a  solution  is  unique  and 
9t-measurable  on  (0,T,) .  Because  tt  is  CADLAG  and  * 
is  measurable,  Tj  is  5  -measurable.  Then,  by  the 
”l 

definition  of  a  Poisson  random  measure  (Jacod  and 


Shiryaev,  1987,  pp.  65-66)  u  is  9  -measurable  =» 

T1  T1 

E„  “  E  +  ,u  )  is  9  -measurable  and,  due  to 

T^  T2-  Tl_  Tl  T1 

A '.3.  E  60  =■  Pathwise  uniqueness  holds  true  on 

*  1 

[O.Tj)  and  60.  Due  to  the  latter,  we  can  repeat 

the  procedure  to  show  that  pathwise  uniqueness  holds 
true  on  rT,,T7)  and  E  60,  and  so  on  for  the 

countable  sequence  of  intervals.  Q.E.D. 


The  interesting  aspect  of  proposition  2.1  is,  that 
the  coefficients  of  (l)'s  fourth  right  hand  term  may 
be  discontinuous  in  E.  This  is  exactly  what  we  need, 
to  construct  a  class  of  hybrid  state  Markov 
processes  that  is  larger  then  the  class  of  solutions 
of  systems  with  Markovian  switching  coefficients. 

The  first  step  towards  this  construction  is 
replacing  *(E,U)  by 

*'<E,u)  -  *(E,U)  X<  [U1<AU)J  U  (F(E)*0)  ),  (2. a) 

where  F  is  a  measurable  mapping  of  Rn  into  (0,1),  * 
and  A  are  measurable  mappings  of  appropriate 
dimensions,  while  the  range  of  A  is  R+.  With  this 
(1)  becomes 

dEt  -  a(Et)dt  +  8(tt)dwt  +  R-iRd  *Ut-'u>  q(dt,du)  + 
+  R-*xRd  P(dt,du).  (2.b) 

Assumptions 

A. "  Define  0'  =  (E6Rn;  F(E)-0), 

(E+#(E,u)j  6  O',  for  all  E6Rn,  u6Rd+1. 

A" ■ 4  Given,  for  all  E6Rn-0'  and  u6R+xRd, 

*<E)-1, 

*<E,u)-*(t,u+Col<l,0,..,0)), 
and  for  any  ken  there  exists  a  constant  M^, 
such  that 

Ti*  “(du)  *  Mk'  for  aU  *EBk- 


A!U5 

3.  A  ( E )  is  on  0'  uniformly  bounded  and  continuous 
in  E . 

b-  (Et)<  t6R+/  exits  0'  at  most  a  countable 
number  of  times. 


jUZ-JIfaggrea 

Given  m(du)  "du^XAi  (du)  and  assumptions  A.  1 .  A.  2  ■  A.  3  . 
A". 4.  A". 5  are  satisfied.  Then  equation  (2.a,b)  has 
for  any  Eg60'  a  pathwise  unique  solution  (Et)- 
Moreover  (Et)  Is  then  a  Markov  process,  of  which  the 
sample  paths  are  measurable  on  the  stochastic  basis 
(0,*,F,P)  . 


Proof : 

Because,  on  O',  A ( E )  is  continuous  in  E  (due  to 
A" ■  5.. a )  and  X(u1<A'),  A'6R,  defines  a  measurable 
mapping  of  R2  into  (0,1)  -  X(u1<A(E))  defines  a 
measurable  mapping  of  RxO'  into  (0,1).  Because  the 
range  of  F  is  (0,1),  we  can  write 

*(  (U!<A(E)]  U  [F(E)*0]  )  -  X(  UX<A(E)  )  V  F ( E ) , 
of  which  both  right  hand  terms  are  measurable.  This 
implies  that  the  supremum  is  measurable,  which 
combined  with  the  measurability  of  *,  makes  that  *' 
is  measurable.  This  ensures  that  (2.b)  is  a  special 
case  of  (1) ,  with  *  replaced  by  *'  according  to 
(2. a).  With  this  we  are  left  to  verify  that  A. 3 . 

A", 4  and  A" . 5  guarantee  that  A' . 3 .  A« .4  and  A' . 5  are 
satisfied,  which  is  straightforward.  Q.E.D. 


Having  theorem  2.2,  we  are  prepared  to  consider  a 
jump  reflecting  boundary  (in  sections  3  and  4)  and 
the  hybrid  state  space  situation  (in  section  5).  But 
first  we  give  a  strong  Markov  characterization  of 
(Et)  If  there  is  no  reflecting  boundary. 


Given  F  vanishes  everywhere  and  the  assumptions  of 
theorem  2.2  are  satisfied.  Then  for  all  t0eRn,  (Et) 
is  a  semi-martingale  strong  Markov  process,  and  its 
extended  generator.  A,  is  given  by: 

At  -  tt  +  1~t  +  ?+f  ,  for  all  f6C2>b(Rn),  (3) 

where 

£f(E)  -i21«i(E)f{i(E)+)(ii2_1[B(E)fl(E)T]ijft  (  (E), 

1  1  (4) 

V^<0)Cf(t+O"f(f,"i£ififE  (C>1 

.  x  <5> 

f  t(()  “Rn£(0)(f(f+O-f(E))  S+(E,dO,  (6) 

and  for  all  Borel  ACRn-(0), 

S'(E,A)  =R-/Rd  *(  *(E,U)6A  )  m(du)  ,  (7) 

S+(E,A)  =  £d  X(  *(E,u)6A  )  dux  x(du).  (8) 

Proof: 

Due  to  A.  3  .  A"  ■  4  ■  A** .  5  and  0'»Rn,  the  9t-predictable 
part  of  Et  ifi 

At  “  £  a  ( Es)ds  +  £  k<'lS~)  £  d  *(Es-,u)  m  (du)  ds . 

Obviously,  (At)  is  of  finite  variation  on  any  finite 
time-interval,  while  {Et~At)  is  a  local  martingale  » 
(Et>  is  a  (special)  semimartingale  (Jacod  and 
Shiryaev,  p.43,  Def.  4.21).  This  immediately  implies 
that  (Et)  is  a  strong  process.  Because  (Et) 

is  a  semimartlngale,  the  generator  A  follows  from 
It6's  differentiation  rule  for  discontinuous 
seeimartingales  (Elliott,  1982).  Q.E.D. 

JP.XZCSSISE.  2EXEBUIEISX1SL  KASKOM,  £BSSESSES 

In  this  section,  we  represent  PD  Markov  processes  as 
solutions  of  an  SDE.  Therefore,  we  consider  (2.a,b) 
with  0-0  and  *  vanishing  on  R~xRd; 

dEt  -  a(Et)dt  VxRd^Ut-'u>- 

,X(  (u^AU/i  U  [F(U*°]  )  p(dt,du)  ,  (9) 

Our  goal  is  to  introduce  a  particular  mapping 
F:Rn-(0,l),  such  that  (9)  has  pathwise  unique 
solutions  which  are  PD  Markov  processes.  The  present 
definition  of  a  PD  Markov  process  (Davis,  1984) 
works  without  such  a  mapping  F.  Instead,  there  is 
given  an  open  subset  0  of  R1*,  with  a  jump  reflecting 
boundary  30,  such  that  (Et)  instantaneously  jumps 
into  the  interior  of  0  just  before  it  would, 
otherwise,  cross  or  travel  through  30.  For  the 
definition  of  a  PD  Markov  process  from  (9)  an 
appropriate  P.  has  to  be  constructed  from  0  and  a. 

The  construction  of  F  will  be  based  on  the  following 
differential  equation,  on  (0,">)xRn, 

dE't  -  a(E't)dt,  t£(0,»),  (10) 

which  has  pathwise  unique  solutions,  assuming  that  a 
satisfies  conditions  and  A. 2 ■  From  this,  we 
define  as  the  set  containing  all  elements  of  30 
that  are  directly  accessible  by  (E't)  from  0: 
in  5  (E630  ;  3  r6(0,«)  and  E ' 060  such  that 

E'r“E  a  E ' t  —6 0 ) .  (11) 

Next  we  introduce  the  following  distance  function, 
<i,(E,i£)  s  inf  (r20  ;  E'0“E  A  E'T6££),  (12) 

which  is,  under  the  above  mentioned  conditions  on  o, 
a  measurable  mapping  of  Rn  into  R.  With  this  we 
define,  for  i6M, 

0 1  -  (E60  ;  da(E,Zfi)  2  1/i),  (13) 

which  are  then  Borel  sets,  and  which  form  the  Borel 
set 

O'* fa  0t.  (14) 

Now  we  define  our  particular  F  as  follows: 

F(E)  -  1  ,  if  E 6Rn-0 ' , 

-  0  ,  else.  (15) 

Due  to  the  above  construction,  F  is  measurable,  by 
which  theorem  2.2  yields: 

_ CgJEllary 

Given  an  open  subset  0  of  Rn,  and  a  mapping  F, 
defined  by  (10)  through  (15).  Then,  under  the 
assumptions  of  theorem  2.2,  equation  (9)  has  for  any 


foeo'  a  pathwise  unique  solution  (  C .  Moreover, 

<  1 1 )  is  then  a  Markov  process,  of  which  the  sample 
paths  are  measurable  on  the  stochastic  basis 
(o,»,r,p) . 

Next,  we  come  to  the  main  result  of  this  section, 
which  implies  that  (it)  is  a  Piecewise  Deterministic 
Markov  process. 

3-2.  -Theorem 

With  probability  one,  the  process  (E^),  o£  corollary 
3.1,  exits  01)3  0  zero  times  on  (0,®). 

SX2&1  ■ 

By  the  definition  of  F,  all  points  of  p  in  R+  become 
active  as  soon  as  (it)  has  e*it  O'.  This  situation 
holds  on  until  ((t)  reenters  O'.  The  reentering  may 
occur  due  to  drift  or  due  to  a  jump  generated  by  a 
point  of  p  in  R+.  Obviously,  the  cases  that  (it) 
--centers  0'  by  drift  without  exit  of  0U££  do  not 
cause  any  difficulties.  In  all  ether  cases,  the 
probability  of  exit  0U3  0  by  drift  is 

|  exp(-s/r|  ds  ■  r  exp(-t/r), 

with  t«inf(l/i  ;  i6N)  and  1/r  the  intensity  of 
points  of  p  in  R+.  Because  (Ef)  exits  O'  at  most  a 
countable  .-.umbel  of  times,  the  probability  of  exit 
0U30  at  least  once  is  then  r/z  exp(-Z/r).  If  all 
points  of  p  in  R+  are  active,  then  because  ten, 

rfS  r,t  exP<"*:/'r)  “  °< 

which  means  a  zero  probability  to  exit  0U£J2.  on 
(0,®).  Q.E.D. 

3  ■  3.  Theorem 

The  process  ( e t ) ,  of  corollary  3.1,  is  a 
semimartingale  strong  Markov  process,  and  its 
extended  generator,  4,  is  given  by: 

At  «  tt  +  ?+f  ,  for  all  £€3>(A) , 
where  t  and  ?  +  are  given  in  proposition  2.3  with 
3=“0 ,  while  the  domain  of  4  is: 

HA)  -  (f  e  c1>b(o)ncb(oui2) ;  7+f(E)«o,  ail  E6£0). 
ExasI: 

Define  a  process  At  as  follows: 

At  -  l  «(es)ds  +  l  x(es_eo')  A(£s_)£d  *(ts-<u)- 

. m(du)ds  +  ^  ^  £d  iMiSi_.u)  du^stdu), 

with  S\  the  ?t-adapted  times  that  (It)  jumps  from 
F.n-0<  into  O',  iil  and  S060, 

Si  “  6^0  >  Si-1  is-6Rn-°'  A  *seo'  >• 

Obviously,  (At)  is  of  finite  variation  on  any  finite 
time-interval,  while  (Et-At)  is  a  loc-l 
?t-martingale.  Subsequently,  (it)  Is  a 
semimartingale.  Application  of  Ita's  differentiation 
rule  for  discontinuous  (piecewise  deterministic) 
semimartingales  to  f(Et)<  with  f  6  C1,  yields: 

f(it>  -  f(t0)  +  l  jfr  f<*s->  Cd«s)i  + 

+  ocist  R+xR<i  -  *<«.->  - 

-  jfr  £(ls-)  (ii(is-'u)i]  P((s),du), 

up  to  indistinguishafeility . 

Substitution  of  A" . 4 . 

p(ds,du)  -  q(ds,du)  +  dsxm(du), 
d(g  -  dAg  +  d(local  martingale), 
m<du)  -  du.XKidii)  , 

and  using  f  e  Cl'b(0)  n  cb(0U£2) ,  yields 

f(it)  -  £(*o)+1E1  l  7~ f(is)(«ces)]ids  +  l  X(ts_60') 

•A(^s')  Id  (f(«s-+*(E  8-.“))  -  f(«s-))  dsxdulXii  (du)  + 

+  ill  l  ld  tf<V+*(V'U))  '  C(tsi-11  + 

+  d(local  martingale), 
up  to  lndistinguishability . 

Next  we  use  the  property  that 

?+f(E)  -  0,  all  Ee££. 

Because  a  is  of  linear  growth  and  (Ct)  is  locally 
bounded,  (a(tt))  is  locally  bounded.  This  implies 
that  (tt)  does  not  increase  while  travelling  through 
0-0'  to  ££,  as  this  takes  a  time  interval  of  zero 
duration.  The  latter  and  the  assumptions  that 
fecb(oui£)  and  f+f(E)-0  for  all  {€££,  imply  that 
f  +  f(Es)-0  for  all  tseo-0'.  With  this, 


f(if)  “  £(*0>  +  ^  ff(EB)  ds  +  d( local  martingale)  + 

+  l  tfa6+^U6,u))  -  f(Es)]  dsxdu.x-fdu). 

Substitution  of  3*  yields 
£(tt)  f(E0)  +  J  4t(Es)  ds  +  d(local  martingale), 

which  implies  that  (Et)  is  a  strong  Markov  process 
with  extended  generator  ( A ,  8(w)].  q.E.D. 

4-. _ PIECEWISE  DIFFUSION  MARKOV  PROCESSES 

Having  obtained  PD  Markov  processes  as  solutions  of 
an  SDE,  the  next  step  is  to  include  diffusion. 
Therefore  we  consider  the  following  SDE: 

dtt  -  a(Et)dt  +  fl(tt)dwt  +R+/Rd  *Ut-'u)  • 

,X(  (u^Mt))  U  [ F ( E ) *0 }  )  p(dt,du),  (161 
which  corresponds  to  (II. a, b)  if  *  vanishes  on  R'xRa 
Initially  we  assume  that  fi(t)8(t)^  is  positive 
definite  for  all  (6Rn,  but  relax  this  assumption 
further  on. 

Now  we  construct  F,  starting  from  the  following 
differential  equation,  on  (0,®)xAn, 
d{'t  -  a(t'»)dt  +  8(t't)dwt,  t£(0,®),  (17) 

which  has  pathwise  unique  solutions  under 
assumptions  A. 1  and  A. 2 .  and  which  defines  a  family 
of  homogeneous  Markov  processes  with  a  measurable 
transition  function 

P't(r,A)  =  P(f 't£A|E'q«E),  all  Borel  A.  (18) 
Because  is  positive  definite,  any  element  of  JO 

is  accessible  by  (E't)  from  0.  Therefore  we 
initially  use  the  following  Euclidean  distance 
function, 

dfl(E,30)  =  inf  ( I E-y|  ;  yeao),  (19') 

which,  obviously,  is  a  measurable  mapping. 

Next,  we  define  the  Borel  sets  0 ^  as  follows, 

0:  -  (EGO  ;  dfl(E,30)  2  1/i),  iGN,  (20) 

and  from  this  the  Borel  set 

°'*i&  0i.  (21) 

As  before,  we  define  our  particular  F  as  follows: 

F(C)  -  1  ,  if  E£Rn-0', 

“  0  ,  else.  (22) 

Obviously,  F  is  measurable,  by  which  theorem  2.2 
yields: 

1U _ Corollary 

Given  an  open  subset  0  of  r.n,  and  a  mapping  t , 
defined  by  (17),  (18),  (19'),  (20),  (21)  and  (22). 

Then,  under  the  assumptions  of  theorem  2.2,  equation 
(16)  has  for  any  Eq6®'  a  pathwise  unique  solution 
{Et>-  Moreover,  (E^.)  is  then  a  Markov  process,  with 
sample  paths  being  measurable  on  the  stochastic 
basis  (c,»,r,P). 

Next,  we  come  to  the  characterization  of  the 
boundary  behaviour  and  the  strong  Markov  property  or 
(«t>- 

Ac2 _ iheacfiB 

With  probability  one,  process  (Et)»  f  corollary 
4.1,  exits  0UJ0  zero  times  on  (0,®). 

Proof: 

By  the  definition  of  F,  all  points  of  p  become 
active  as  soon  as  (E^)  has  exit  O', say  at  moment  T, 
which  situation  continues  until  (E^l  has  reentered 
O',  say  at  moment  T+d.  The  exit  may  occur  due  to 
diffusion  or  due  to  a  jump  generated  by  a  point  of  p 
in  R+.  Obviously,  the  cases  that  (Ep)  exits  0-0'  by 
diffusion  without  entering  30  do  not  cause  any 
difficulties.  In  all  other  cases  we  know  from  the 
proof  of  theorem  3.3  that  d  has  an  exponential 
distribution  of  which  both  the  mean  and  the  standard 
deviation  equals  r-0+.  with  this,  it  follows  that, 
for  any  E€0',  the  probability  of  entering  and 
exiting  30  within  1/t  is: 

r-1  P'j (r,Rn-0-30)  £  r*1  PE'(r,(ySRn  ;  |E~y|  >  <  )) 
with  c-inf  { l/i  ;  ieN)  . 

Because  (E^)  is  a  diffusion  and  t>0,  the  right  hand 
side  is  of  order  r  (Gihaan  and  Skorohod,  1972,  p. 

64).  As  this  situation  may  occur  a  countable  number 
of  times,  we  have  to  divide  by  r,  yielding  order 
(r/t),  of  which  the  limit,  rio,  is  zero.  Q.E.D. 

lxJ _ Proposition 

Given  the  assumptions  of  theorem  4.2  are  satisfied. 
Then  for  all  f0€0',  { E t )  is  a  semimartingale  strong 
Markov  process,  and  its  extended  generator,  A,  is 


given  by:  - 

At  -  it  +  }+t  ,  for  all  tev(A), 
where  t  and  are  given  in  proposition  2.3,  while 
the  domain  of  A  is: 

v(a)  -  (f  e  c‘<5(0)ncb(ouao) ;  ?+f(E)-o,  all  Eejo). 

Proof :  Similar  to  the  proof  of  proposition  3.3, 
except  that  now  }+t(ts)-Q,  tor  all  Es£0-0',  follows 
from  fec(0U30) .  Q.E.D. 

Finally,  we  consider  the  more  general  situation  with 
8(E)8(E)T  being  positive  semidefinite.  The 
construction  of  F  works  according  to  equations  (17), 
(18),  (20),  (21)  and  (22),  but  with  distance 
function: 

dfl(t,iLC)  -  inf  (r20;  (IS  fl  E,  r)x()  ),  (19) 

where  IS  Is  the  subset  of  JO  tnAt  is  accessible  by 
(t't)  from  0,  ( )  is  the  empty  set  and  E,  _  is  the 
closure  of  an  n-dimensional  ellipsoid,  with  centre 
E+«(E)r  end  shape  defined  by  covariance  fl(E)8(E)*r. 
Obviously,  d g(.,iS)  is  measurable,  by  which  the  Oj/s 
and  O'  are  Borel  sets  and  F  is  measurable,  and  we 
get : 

_ Corollary 

Given  an  open  subset  0  of  Rn,  and  a  mapping  F, 
defined  by  (17)  through  (22).  Then,  under  the 
assumptions  of  theorem  2.2,  creation  (IV)  has  for 
-ny  Eq60'  a  pathwise  unique  solution  {1^).  Moreover, 
lip)  then  a  Markov  process,  with  sample  paths 
being  measurable  on  the  stochastic  basis  (fl,?,F,P). 

Next,  we  come  to  the  main  result  of  this  section. 

4 ■ 5  Theorem 

with  probability  one,  the  process  (E^),  of  corollary 
4.4,  exits  0U£S  zero  times  on  (0,®). 

Proof : 

By  the  definition  of  F,  all  points  of  p  in  R+  become 
active  as  soon  as  (Et)  has  exit  O'.  This  situation 
holds  on  until  (Er)  reenters  O'.  The  reentering  may 
occur  due  to  drift  and/or  diffusion  or  due  to  a  jump 
generated  by  a  point  of  p  in  R+.  obviously,  the 
cases  that  (Et)  reenters  0'  by  drift  and/or 
diffusion  without  exit  of  OUJJ)  do  not  cause  any 
difficulties.  Of  those  cases  where  £S  Is  accessible 
through  drift  only,  we  follow  the  proof  of  theorem 
3.1.  Say  aoa  is  the  subset  of  £S  that  can  only  be 
entered  by  (E't)  due  bo  drift.  For  all  other  cases 
we  then  notice  that  a  strictly  positive  type  (19) 
distance  d.  at  the  moment  of  exit  O',  corresponds 
with  a  strictly  positive  Euclidean  distance  from 
3  0-3  0„ ■  due  to  the  local  boundedness  of  |a(it)|  and 
|S(lt)|.  Subsequently,  we  may  follow  the  proof  of 
theorem  4.2  for  these  cases.  Q.E.D. 

1.6  .-Tfaggcerc 

Given  the  assumptions  of  corollary  4.4  are 
satisfied.  Then  for  all  EnEO',  (Et)  is  a 
semimartingale  strong  Markov  process,  and  its 
extended  generator.  A,  is  given  by: 

At  -  tl  +  7+f  ,  for  all  f£D (A)  , 
where  t  and  are  thos»  given  in  proposition  2.3, 
while  the  domain  of  A  is: 

V(A)  -  (f  6  CJ'b(0)ncfa(0U££) ;  7  +  f(E)-0  all  *€££). 

Proof :  Similar  to  the  proofs  of  theorem  3.3  and 
proposition  4.3. 

5^  THE  HYBRID  STATE  SPACE  SITUATION 

In  this  section  we  explicitly  consider  the  hybrid 
state  space  situation  for  a  system  of  the  form 
( 2 . a , b) ,  in  such  a  way  that  there  is  no  need  of 
assuming  a  particular  i  or  A.  As  such,  all  jump 
reflecting  boundary  results  of  the  former  sections 
fit  into  the  results  of  this  section.  For  ease  of 
notation  and  interpretation,  we  rewrite  the  SDE  form 
(2.a,b)  by  replacing  the  Poisson  random  measure,  p, 
by  a  multivariate  counting  process,  »t,  such  that 
the  pathwise  uniqueness  of  (2)'s  solution  is 
preserved.  We  do  that  by  defining,  for  all  Borel 
UCR+xRd, 

'tfu)  ~l  {,  *<•  (U1<*(ES_)  J  U  [F(Es_)*0J  )  p(ds,du), 

and  then  rewriting  (II)  as  ' 

dtt  ■  a  f  E  t) dt  +  S(tt)dwt  +R-^Rd  *(Et-,u)  <?(dt,du)  + 


+  R+xJtd  d«t(du)  •  (23. b) 

The  main  objective  of  this  section  is  to  show  that 
the  last  term  of  (23. b)  generates  a  particular  type 
of  jump:  a  jump  in  (it)  that  anticipates  a 
simultaneous  switching  of  U1*).  For  short  we  refer 
to  this  type  of  jumps  as  hybrid  lumps.  Notice  that 
these  hybrid  lumps  are  in  some  sense  unexpected,  as 
all  coefficients  of  (23. a, b)  are  non-anticipating. 

To  show  these  hybrid  jumps  explicitly,  we  need  some 
preparation. 

5^-1 _ ksama 

Under  assumptions  A . 1 .  A. 2.  A". 3.  A" .4  and  A" ■ 5.  the 
pair  of  equations  (23. a, b)  has  for  any  E0e0'  a 
pathwise  unique  solution  (Et'l't)'  wbere  jit  a 
multivariate  counting  process  on  R+xR+xRd  of  a 
predictable  intensity,  At=A(Et_).  Moreover  both 
(*t'yt)  and  l S t >  are  then  semrmartingale  strong 
Harkov  processes,  of  which  (Et)  is  indistinguishable 
from  the  one  in  theorem  2.2. 


Proof: 

It  follows  from  theorem  2.2,  that  the  system  of 
equations  (2.a,b)  and  (23. a)  has,  for  any  Borel  U,  a 
pathwise  unique  solution  <  E t • r  t  <u>  > .  With  this, 
system  (2.a,b),  (23. a)  has  a  pathwise  unique 
solution  (Et.'t)-  Obviously  all  potentially  active 
points  of  p,  that  are  in  R+xR+xR®,  are  collected  by 
>t  (n  a  predictable  way,  by  which  we  can  write 
R+xRd  *'(Et-'u>  *<  (Ui<AO:t_)l  U  (F(Et_)*oj  ). 

-P(dt,du)  -  R+/Rd  *(Et_,u)  dvt(du)  , 

up  to  indistinguishability.  This  implies  that  the 
solution  of  (2-b)  is  indistinguishable  from  the 
solution  of  (23. b).  Q.E.D. 


Now  we  are  prepared  to  consider  the  hybrid  state 
space  situation.  Therefore  we  assume  that  the  first 
component  of  E t  is  H-valued,  with  HCN= (1,2,..),  and 
that  we  can  write  the  first  scalar  equation  of 
(23. b)  as  follows: 

d^t  “  jt'jl  dvt(du),  (23. c) 

with  a  mapping  of  RnxR+xRd  into  the  integer 
lattice,  Z. 

Next  we  assume  that  *  satisfies,  for  all  A (( )  ] , 


F(E,U)%^  X 


Ho x(i'E)  s  ui  c  iIox(i'E)]  * 


(1,E,U) 


where  »  is  a  measurable  mapping  of  KxRnxRd  into 
ZxRn~l,  and  1  is  a  measurable  napping  of  NxRn  into 
R+,  such  that  l(i,.)-0  for  all  i£N-M,  and 

iiN  A<<>  • 

Mcjover,  we  assume  that  for  all  n£N,  E6ZxRn-l  and 
U6Rd, 

*l(1,E,U)  “  1-E  i  ,  (25) 

which,  together  with  (24)  and  x(i,.)~0  for  all 
i£N-M,  implies  that  if  E  q£M,  then  ( E 1 t )  in  R+xM. 
Substitution  of  (24)  and  (25)  in  (23.b,c)  and 
subsequent  evaluation  yield 


df1! 


t  -  i*  niK  x  HoX(i'ft-)  *  u*t.  <  J0X(i,Et_) 


2. 


.  .(4-E1t-)  d't(dUlxRd) ,  (26. a) 

with:  u  s  =  u1-kA(Es), 

for  some  integer  k  such  that  0  <  u*6  £  a(Ee), 
dIt  “  fl(tt>dt  +  £(Et)dvt  tg-J^d  t(«t-'u)  q(dt,du)  + 

+  £d  *U1t'tt-'U)  d»t(R+xdu),  (26. b) 


't(W)  “  ^  K  (u1£A(Es.)]U(F(Es.)"0)  )  p(ds,du). 


1  H  (26. C) 

all  Borel  UCR^xR0,  where  underlining  of  a  vector 
refers  to  all,  but  the  first,  components  of  that 
vector. 


Assumptions 

A  ■  *■  Given,  for  all  EGRn-0'  and  u6R+xRd, 

ME)  -  ,^(M)  -  1. 

*i(1,C,U)  -  1-El, 
m(du)  -  dujxxfdu). 

For  all  k£N  there  exists  a  constant  HR,  such 
that  for  all  teBu. 

1(1, E)  (11-fil  *  £d  I*(1,E,U)I  a  (dU)  ]  £  M*. 


a-  For  all  E60' , 

l(i,E)  is  continuous  in  E, 
i (i,E)“0  for  all  ieN-H, 

A(E)  =  1(1,1)  is  uniformly  bounded. 

b-  (£(■),  t£R+,  exits  0'  at  most  a  countable  number 
of  times. 

iJi — Theorem 

Given  the  hybrid  space  O'  s  0'n(MxRn_1). 

Under  assumptions  AjJ.  through  Ai£,  the  system  of 
equations  (26.a,b,c)  has  for  any  EqEO'  a  pathwise 
unique  solution  (Ct,»^).  Moreover  (Et)  is  then  a 
semimartingale  strong  Markov  process  in  R^O'. 

Proof: 

Due  to  A .  3  and  A ■  5 ■  a  .  (24)  defines  Ms  a  measurable 

mapping  (see  proof  of  theorem  2.2),  by  which 
( 26 . a , b, c)  is  a  special  case  of  (23.a,b,c).  Next  we 
show  that  Aj^i  implies  A" .  4  .  by  which  lemma  5.1  and 
(24)  imply  that  the  solution  of  (23.a,b,c)  is 
indistinguishable  from  the  solution  of  (26.a,b,c). 

To  arrive  at  A" ■ 4 .  we  start  from  A ■ 4  and 
subsequently  use  A . 5 . a .  interchange  order  of 
integration  and  substitute  (24).  Q.E.D. 


Due  to  its  extensive  form,  equation  (26.a,b,c)  hides 
the  results  for  which  the  above  analysis  has  been 
carried  out.  Therefore,  we  take  a  closer  look  at  it 
in  case  that  p  has  no  points  in  R“.  Then,  (26. b) 
becomes 

dlt  -  E(lt)dt  +  2<(t)dwt  +  Id  2LU1t'tt-<U)d''t(R','*dU) 

(27. a) 

Moreover,  to  avoid  the  use  of  equations  (26. a, c),  we 
go  over  to  the  common  descriptive  way  of  formulating 
(vt)  and  (  E  1t  > : 

{  v  j_ )  is  a  multivariate  counting  process 
characterized  by  the  9t-predictable  intensity,  Tt, 
rt  '  *<lt->  t1  +  F<*t->  Ifg  1/T),  (27. b) 

and  a  deterministic  jump  measure  w(du). 

(Mt)  is  a  process  with  a  countable  state  space,  N, 
and  with  an  9t-predictable  rate,  r^  t,  of  jumping 
from  lit_=j  to  Mt=i,  i*j,  J' 

rij,t  =  Mi,(j,it-I)  (1+  F(.,It_)  l|g  1/r],  (27. c) 
while  rlj(t  S  rt. 

From  this  formulation,  we  easily  notice  the 
interesting  effect  that  E^t  appears  in  the 
coefficient,  v,  of  (27. a) 's  third  right  hand  term. 
This  means  that  ®  ( E  V  ,  E  t- <  U)  anticipates  a  switching 
from  Mr-  to  '  t'  and  thus  a  jump  of  (It) 
anticipates  a  simultaneous  transition  of  (E^t). 
verify  that  the  anticipating  coefficient  £  already 
appears  in  (26.a,b,c),  while  there  is  no 
anticipating  coefficient  in  equation  (23.a,b,c).  As 
the  solutions  of  both  equations  are 
indistinguishable,  we  conclude  that  (23.a,b,c)  is 
the  canonical  representation  of  a  system  with  hybrid 
lumps .  while  (26.a,b,c),  with  the  anticipating 
coefficient,  is  the  representation  that  is  more 
useful  when  it  comes  to  the  realization  of  Markov 
models  with  hybrid  lumps. 

Remark:  If  X  ( • ,  (( j,,!) )  is  t-invariant,  then  (E^t)  *-s 
a  countable  state  Markov  process.  In  this  case 
(27. a)  can  straightforwardly  be  obtained  from  a 
classical  system  like  (1)  of  which  all  coefficients 
are  continuous.  For  the  situation  that  (It) 
continuous,  i.e.  v-0,  see  Brockett  and  Blankenship 
(1977).  For  some  applications  with  hybrid  lumps, 
i.e.  **0,  see  sworder  (1972),  Blom  (1984)  and 
Mariton  (1987). 
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